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I.  GENERAL  INTRODUCTION 


The  Marangoni  phenomenon  refers  to  the  fluid  motions  induced  in  the  vicinity  of  a  surface 
of  separation  between  two  different  fluids.  At  this  surface,  the  thermodynamic  properties 
(such  as  density,  pressure,  composition,  ...)  undergo  rapid  variations  on  a  very  small  scale 
(the  interfacial  phase  thickness).  The  experimentation  on  such  interfaces  in  equilibrium 
configurations  shows  that  the  result  of  these  rapid  variations  can  be  characterised  by  an 
internal  energy  per  unit  area,  that  is  commonly  referred  to  interfacial  tension  (for  liquid/gas 
interfaces,  the  name  surface  tension  is  generally  preferred).  This  allows  to  consider  the 
fluid/fluid  interfaces  as  discontinuities  in  the  above-mentioned  variables.  Since  the  interfacial 
tension  is  in  turn  a  function  of  intensive  thermodynamic  variables,  interfacial  tension 
variations  along  the  surface  can  be  present,  provided  that  gradients  of 
temperature/composition  exist.  These  surface  tension  gradients  (tangential  stresses)  may  then 
create  motions  at  the  interface,  that,  due  to  viscosity,  generally  extend  relatively  for  into  the 
bulk  of  the  adjacent  liquids. 

This  phenomenon,  also  called  thermocapillary  convection  (when  the  driving  gradient  is 
thermal),  or  solutocapillary  convection  (when  it  is  compositional)  may  result  in  extremely 
different  forms  of  convective  motions,  such  as  stationary  motions,  oscillations  (waves),  or 
even  turbulent  behaviours  (interfacial  turbulence).  These  different  kinds  of  motions  are 
generally  divided  into  two  broad  classes,  according  to  the  direction  of  the  gradients  with 
respect  to  the  interface. 

The  first  case,  called  Marangoni  convection,  occurs  when  the  gradients  are  parallel  to 
interface.  Since  a  stress  is  created  along  the  interface,  at  that  no  counteracting  force  is 
present  when  the  fluid  is  at  rest,  motions  set  in  whatever  the  intensity  of  the  gradient.  A 
steady  state  is  generally  reached  when  the  hydrodynamic  pressure  gradients  created  by  fluid 
motions  are  able  to  balance  the  imposed  surface  tension  stresses.  An  example  of  this 
phenomenon  occurs  when  a  cylindrical  bridge  of  liquid  is  confined  between  two  rigid  disks 
maintained  at  different  temperatures  (this  is  the  half-zone  configuration  often  used  in  the 
studies  of  other  more  complicated  techniques  directed  to  the  processing  of  crystals).  A 
convective  axisymmetric  toroidal  motion  is  generally  observed  (the  motion  is  from  hot  to 
cold  at  the  surface,  provided  that  the  surface  tension  is  decreasing  with  temperature).  For 
higher  temperature  gradient,  other  kinds  of  more  complicated  (non-axisymmetric)  oscillating 
motions  can  be  obtained,  and  are  responsible  for  poor  properties  (surface  state,  composition 
homogeneity,  ...)  of  crystals  grown  by  these  techniques.  Other  technological  processes 
stimulating  the  need  for  an  understanding  of  Marangoni  motions  occur  in  the  drying 
technology  (where  Marangoni  motions  may  at  the  contrary  have  a  positive  role),  in  the 
coating  industry,  in  the  laser  processing  of  compact  disks,  in  heat  pipes  technology,  in 
boiling,  ... 

A  second  situation,  called  Marangoni-Benard  convection,  occurs  when  the  driving  gradients 
are  perpendicular,  rather  than  parallel,  to  the  interface.  In  this  case,  no  tangential  stress  is 
created  (because  the  interface  is  at  constant  temperature/  composition),  and  the  mechanism 
responsible  for  convection  is  intimately  associated  with  the  presence  of  fluctuations 
(thermodynamic,  or  of  external  origin)  of  the  intensive  properties,  coupled  with  the 
hydrodynamic  effects  (transport  of  heat,  or  of  mass).  For  example,  consider  a  fluid  layer 
lying  on  a  heated  rigid  plate,  and  open  to  air.  When  a  small  fluctuation  of  temperature  (say 
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an  increase)  occurs  at  a  point  on  the  free  surface  (it  may  also  originate  in  the  liquid  and 
diffuse  to  the  surface),  the  surface  tension  ^ 

is  slightly  lower  at  that  point.  This  creates 
surface  tension  gradients  along  the  surface, 
producing  motions  driving  fluid  away  from 
the  point.  Continuity  of  the  fluid  then 
requires  fluid  coming  from  the  bulk  phase 
below  the  point,  to  the  interface.  Now, 
since  a  gradient  of  temperature  is  imposed 
by  the  heating,  the  fluid  arriving  at  the 
surface  is  slightly  hotter,  and  this  lowers 
the  surface  tension  again,  thus  amplifying 
the  motions.  This  may  also  result  in  j 
stationary  convective  structures,  that  create  | 
a  tessellation  of  the  layer  in  an  almost 
periodic  way  (see  fig.  1).  An  important  ^ 
difference  exists  between  Marangoni  [ 
convection  and  Marangoni-Benard 
convection.  For  the  former,  motions  exist 
whatever  low  is  the  value  of  the  gradient 
imposed  along  the  interface.  At  the 
contrary,  for  the  latter,  a  critical  value  of  pjgm.g  j  j.  Hexagonal  convective  structure 
the  imposed  gradient  must  be  exceeded  by  heating  a  thin  liquid  layer  heated 

before  motions  occurs.  This  is  associated  below.  The  fluid  moves  upwards  at  the 

with  the  presence  of  both  thermal  gf  ^  hexagon,  and  downwards  at  its 

diffusivity  and  viscosity,  the  first  of  which  pgj.jpheral. 
damps  the  thermal  fluctuations,  and  the 
second  of  which  damps  the  fluid  motions 

created  by  surface  tension  fluctuations.  The  competition  between  both  these  stabilising  effects 
and  the  destabilising  effect  of  surface  tension  variations  is  typically  expressed  through  the 
dimensionless  Marangoni  number 

-(jT-ATh 


where 

Oj  surface  tension  variation  with  temperature 

h  :  thickness  of  the  layer 

/i  :  dynamic  viscosity 

K  :  thermal  diffusivity 

AT  ;  temperature  difference  between  the  rigid  conducting  plate  and  the  free 

surface  in  the  absence  of  convection  (AT=^h,  where  ^  is  the 
imposed  gradient). 

A  theoretical  linear  stability  analysis,  such  as  the  one  presented  in  this  report,  leads  to  the 
result  that  convective  motions  develop  in  the  fluid  provided  that  Ma>Ma^—80  (when  the 
free  surface  Biot  number  is  zero,  meaning  that  the  heat  flux  crossing  it  stays  constant).  The 
linear  stability  analysis  also  provides  the  horizontal  length  scale  of  the  periodic  strucmres 
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(size  of  the  hexagonal  cells  at  the  threshold)  as  approximately  equal  to  3.14  h. 

Thus,  the  fluid  behaviour  undergoes  a  drastic  qualitative  change  when  the  Marangoni  number 
exceeds  the  critical  value,  a  phenomenon  which  is  called  a  bifurcation,  in  the  frame  of 
nonlinear  stability  theories.  The  theoretical  interest  of  the  Marangoni-Benard  instability  relies 
in  the  fact  that  it  provides  a  relatively  simple  and  well-controlled  example  of  an 
hydrodynamic  instability,  that  can  be  studied  both  theoretically  and  experimentally.  Note  that 
in  some  situations,  it  is  very  similar  (on  the  point  of  view  of  theoretical  techmques  used  near 
the  instability  threshold)  to  the  more  classical  Rayleigh-Benard  convection  (that  occurs  in  Ae 
same  experimental  design,  but  when  the  destabilising  effect  is  the  density  variation  with 
temperature,  i.e.  the  buoyancy  force).  However,  important  differences  exist  m  other 
conditions  (f.ex.  when  the  imposed  gradient  is  increased,  and  that  secondary  instabilities  and 
transitions  to  more  complex  time-dependent  simations  occur),  which  justify  the  analysis  of 
Marangoni-Benard  convection  as  an  hydrodynamic  instability  capable  of  displaying 
phenomena  such  as  pattern  formation,  defect  generation  (see  the  distorted  hexagons, 
pentagon-heptagon  pairs, ...  of  fig.  1)  and  motion,  wavelength  adjustment  mechamsms,  effects 
of  lateral  walls  (and  other  experimental  imperfections),  influence  of  symmetries,  resonance 
effects,  transient  phenomena,  secondary  instabilities  and  transitions  to  mrbulence. 

An  important  practical  feature  of  the  Benard  instabilities  is  the  drastic  enhancement  of  heat 
transfer  which  is  created  by  convection.  This  is  generally  measured  by  the  Nusselt  number, 
representing  the  dimensionless  ratio  of  the  total  (convective  and  conductive)  heat  flux  to  the 
conductive  heat  flux  :  its  value  is  thus  always  superior  or  equal  to  one,  and  is  expected  to 
increase  with  the  distance  to  the  threshold  (generally  measured  by  e=(Ma-MaJ/Ma,).  The 
dependency  of  the  Nusselt  number  on  the  heating  conditions  (thus  on  e)  is  not  well-known 
for  high  values  of  the  heat  flux,  and  its  knowledge  is  of  fundamental  importance  in  all 
domains  where  thermal  energy  has  to  be  transported  from  one  place  to  another  through  liquid 
phases.  In  this  respect,  it  has  been  recognised  that  Marangoni  convection  might  be  a 
determinant  factor  in  the  amount  of  boiling  heat  transfer  under  reduced  gravity  conditions. 
The  understanding  of  this  heat  transfer  mechanism  is  also  of  tremendous  importance  to 
optimise  evaporators  which  are  used  in  large  satellites,  as  capillary  pumps  (heat  pipes)  to 
transfer  thermal  energy  from  the  place  where  it  is  dissipated  (electronic  boxes)  towards 
radiators. 

In  some  problems,  the  Nusselt  number  may  have  to  be  replaced  by  other  quantities,  which 
provide  a  better  measure  of  the  enhancement  of  heat  transfer  created  by  convection.  The  heat 
flux  may  for  example  be  fixed  as  the  external  parameter,  and  we  will  be  interested  in 
calculating  the  cooling  of  the  heat  dissipating  component  by  convection.  This  latter  approach 
will  be  adopted  in  this  report,  as  explained  in  section  III,  where  a  nonlinear  analysis  of  the 
convection  problem  is  presented,  leading  to  estimates  of  the  so-called  bulk  temperature 
decrease.  Other  results  are  obtained  concerning  the  important  problem  of  the  selection  of  the 
wavelength  of  the  steady  convective  structure  at  given  Marangoni  number.  In  section  II,  we 
begin  by  defining  the  important  characteristics  of  the  Marangoni-Btad  instability  (the 
dimensionless  parameters  entering  the  problem),  and  investigate  about  the  modelling  of  the 
effect  of  evaporation  and  surface  deformation.  Some  classical  results  will  also  be  recalled, 
and  the  limits  of  applications  of  our  model  will  be  determined.  Finally,  numerical  simulations 
(Galerkin  method)  and  comparison  with  analytical  results  of  section  III  will  be  presented  in 
section  IV. 
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11.  STATEMENT  OF  THE  PROBLEM  -  LINEAR  STABBLITY  ANALYSIS 

11. 1.  Rate  of  evaporation  in  the  presence  of  surface  deformation 

Consider  a  layer  of  a  pure  liquid,  lying  on  a  heated  rigid  plate  (z=0)  and  in  contact  with  its 
own  vapour  at  a  free  interface  whose  shape  is  given  hy  z=h(x,t),  in  which  x  is  the  horizontal 
coordinate  and  t  is  time.  At  this  interface,  evaporation  may  occur,  such  that  a  mass  flux 

J  =  7.7fis  allowed  (n  =  {-h' ,\)IN  is  the  free  surface  normal  pointing  to  the  vapour, 
is  the  normalisation  factor,  and  a  prime  denotes  a  derivative  with  respect  to 
x).  In  many  smdies  of  evaporating  liquid  layers,  the  mass  flux  J  is  assumed  to  obey  the  well- 
known  Hertz-Knudsen  equation,  derived  from  the  kinetic  theory  of  perfect  gases.  This  is 
valid  for  a  flat  interface  along  which  the  equality  of  liquid  and  vapour  temperatures  is 
assumed,  and  reads 

/  =  /3 

where  /3  is  the  constant  accommodation  coefficient,  M  is  the  molecular  weight  of  vapour, 
p^(r)  is  the  saturation  pressure  at  surface  ternperamre  T,  Po(r)  is  the  vapour  pressure  just 
beyond  the  interface,  and  R  is  the  universal  gas  constant. 

In  this  work,  in  order  to  start  from  a  more  general  form  of  the  nonequilibrium  mass  flux 
equation  which  should  in  particular  be  valid  for  a  moving  interface  of  arbitrary  shape  (and 
for  any  form  of  the  vapour  state  equation),  we  suggest  to  use  the  macroscopic  approach,  as 
was  proposed  in  [1].  First,  note  that  the  Clausius-Clapeyron  equation  can  be  derived  from 
the  equality  condition  iJL,(p,,TJ=fii(PiJ)  between  vapour  and  liquid  chemical  potentials. 
Developing  this  formula  by  the  use  of  classical  thermodynamics  relations,  the  Clausius- 
Clapeyron  equation  gives  the  slope  of  the  coexistence  curve  : 

dp  _  PvP^  (2.2) 

If  iPrP)T 

obtained  assuming  thermal  equilibrium  Ty=Ti(=T),  and  mechanical  equilibrium ;7y=/7,(=pj. 
In  this  relation,  and  Pi  are  respectively  the  vapour  and  the  liquid  volumic  masses,  and  L 
is  the  latent  heat  of  evaporation.  Since  this  mechanical  equilibrium  condition  is  not  valid  in 
general,  the  equality  of  chemical  potentials  does  not  hold  anymore.  The  difference  of 
chemical  potentials  is  a  generalised  thermodynamic  force  giving  rise  to  a  nonequilibrium 
mass  flux  across  the  interface.  This  obviously  suggests  the  use  of  the  thermodynamics  of 
irreversible  processes  [2],  through  the  phenomenological  law  : 

J  =  K[{jLi{pi,T)  -  ii^(Py,T)]  (2.3) 

where  K  is  the  positive  phenomenological  coefficient.  Note  that  the  hydrodynamic  definition 

of  the  mass  flux  J  is  _ _ _ 

J  =  P,(v,-v^.n  =  p,(v,-Vj.).«  (2-4) 

where  "v; ,  and  are  respectively  the  vapour,  liquid  and  interface  velocities. 

Note  that  if  the  ternperamre  jump  is  not  neglected,  then  other  phenomenological  coefficients 
appear,  since  there  exist  also  a  "thermal  force"  AT=TrT,.  We  discard  this  case  in  the 
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presented  analysis. 


Let  p,(r)  be  a  function  defining  the  equilibrium  (saturation)  pressure  at  a  given  temperature 
T  (for  example  the  Clausius-Clapeyron  coexistence  curve,  or  any  fitting  of  experimental 
points).  Assume  that  the  liquid  and  vapour  state  equations  are  respectively 

p,  =  p?(^,D  ,  p,  =  pl(p,T) 


and  that  the  inequalities 

\prpm  ^ ,  \p.-pm  ^ , 
p,a)  ’  p,a) 

are  valid  everywhere  along  the  liquid-vapour  interface.  Then,  using  Taylor  expansion  around 
the  point  (p,(T),T),  at  which  chemical  potentials  are  equal  (by  definition  of  p,(T)),  we  get  a 
linearised  form  of  eq.  (3)  as 


J  =  K(T) 


Py{T)  pfj) 

where  p.fT),  p,(7)  and  K(r)  stand  for  p>,(T),r),  p^(p,(r),T) 


Ps<T>-Pv  Ps(^-Pi 


(2.5) 

and  K(p,(T),T)  respectively. 


Note  that  if  p^-pi  in  eq.  (5),  we  get  a  form  J=K  (p,(D-p,)/py(r),  presenting  a  formal 
analogy  with  the  Hertz-Knudsen  equation  (although  obtained  independently),  and  leading  to 
a  rough  estimate  of  the  phenomenological  coefficient  ^:=|8pv(r)  (M/licRTf^  Of  course,  only 
experimental  measurements  of  K  can  lead  to  satisfactory  values. 


However,  the  approximation  Pv~~Pt  is  only  valid  for  weakly  curved  interfaces,  and  when 
dynamical  effects  can  be  neglected  in  the  momentum  balance  at  the  interface.  A  more  general 
form  of  equation  (5)  can  be  obtained  by  considering  the  Laplace  pressure  difference 
p-p^=aH  (H  is  the  mean  curvature  of  the  interface)  in  eq.  (5),  thus  leading  to  the  equation 


J  = 


pfTf-p^T) 

Pi(T)pJJ) 


K(X) 


P,{T)-p^-a 


pSX) 


P,(T)-p^(T) 


H 


incorporating  the  effect  of  surface  tension  on  the  rate  of  phase  change. 


(2.6) 


Note  that  eq.  (6)  is  only  valid  provided  \p,-p,(T)\  <  \aH\  <p,(T).  In  pure  weightlessness, 
every  isothermal  liquid-vapour  system  with  an  interface  of  constant  curvamre  H=Ho  is  in  a 
state  of  mechanical  equilibrium.  Since  no  mass  flux  occurs  across  the  interface,  eq.  6 
converts  into 


Pv 


=  Ps(6)  -  a 


P,(0)-Pv(0)  ' 


where  0  is  the  constant  temperamre  of  the  isothermal  system,  this  equation  defines  a 
connection  between  the  samration  vapour  pressure  Ps(0)  at  a  plane  phase  interface  and  the 
vapour  pressure  py  at  an  equilibrium  interface  of  curvamre  Hg.  In  the  case  of  spherical 
interfaces,  this  is  known  as  the  Thomson  equation. 
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II.2.  Calculation  of  the  rate  of  evaporation  in  the  basic  state 

In  this  section,  we  compute  the  rate  of  evaporation  of  a  liquid  layer  of  depth  h  lying  on  a 
heated  rigid  plate  maintained  at  a  constant  temperature  T5.  The  system  is  assumed  as 
unidimensionnal  (all  the  variables  only  depend  on  the  vertical  coordinate  z,  and  the  only  non¬ 
zero  component  of  the  velocity  in  both  liquid  and  gas  phases  is  the  vertical  one,  that  will  be 
noted  w).  The  vapour  phase  is  infinitely  deep.  The  solution  of  the  hydrodynamic  equations 
that  will  be  obtained  will  be  called  the  basic  solution  (or  the  reference  solution),  and  its 
stability  against  hydrodynamic  two-dimensional  fluctuations  will  be  studied  in  section  II. 3 
(Marangoni  problem). 

In  the  presence  of  evaporation,  the  relation  (4)  indicates  that  if  the  liquid  is  at  rest  (W|=0), 
the  interface  moves  with  a  vertical  velocity  Wj;=-J/p,.  The  second  relation  (4)  then  indicates 
that  the  vertical  velocity  of  the  vapour  is  Wy=J  (l/pv-l/pi)=J/pv-  This  means  that  the  liquid 
depth  will  decrease,  and  that  it  may  finally  completely  evaporate.  In  order  to  obtain  a  steady 
reference  solution  (and  a  steady  mass  flux  J),  we  will  assume  in  the  following  that  fresh 
liquid  is  injected  at  the  bottom  plate  z=0,  at  a  rate  that  exactly  matches  the  rate  of  liquid 
which  is  evaporated  at  the  free  surface  z=h.  This  is  attempted  in  order  to  mimic  the  steady 
regime  solution  which  is  obtained  in  heat  exchanging  devices  such  as  heat  pipes,  where  liquid 
condensed  at  the  condensor  is  continuously  brought  back  to  the  evaporator,  such  that  a  steady 
regime  may  be  obtained.  The  rigid  bottom  boundary  may  thus  be  considered  as  a  porous 
material,  through  which  a  velocity  Wb=J/pi  is  injected.  Now,  the  whole  liquid  layer  thus 
possesses  the  velocity  w,=Wb=J/pi  (because  the  fluid  is  incompressible,  such  that  dw/dz=0), 
and  the  relation  (4)  indicates  that  Wj;=0.  The  interface  is  thus  fixed  in  space  at  the  location 
z=h. 

The  effective  rate  of  evaporation  J  is  given  by  the  relation  (5),  which  is  seen  to  depend  on 
the  vapour  and  liquid  pressures  py  and  pi  on  each  side  of  the  interface,  and  on  the  interfacial 
temperature  T.  Their  value  may  only  be  obtained  by  solving  the  hydrodynamic  equations  in 
both  liquid  and  gas  phases,  with  suitable  boundary  conditions. 

In  both  phases  (indices  v  and  1  are  omitted  for  simplicity),  the  full  Navier-Stokes 
(conservation  of  momenmm)  equations  for  Newtonian  fluids  read 

p(—+  ( V  V)  V)  =  pAv-  Vp  +  pg  (2.7) 

dt 

where  p  is  the  dynamic  viscosity  (assumed  to  be  independent  of  temperature),  andg=  -g  1^ 

is  the  gravity  vector  (T^  is  the  unit  vector  along  the  vertical  z-direction).  As  the  system  is 
unidimensional  and  the  velocities  are  independent  of  z  and  t  (steady  problem)  in  both  phases 
(the  gas  is  also  assumed  to  be  incompressible,  which  is  realistic  for  the  range  of  velocities 
involved),  the  equations  (7)  reduce  to  an  equation  for  the  hydrodynamic  pressure 

±  (2.8) 

dz 
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Now,  we  have  to  consider  the  equation  governing  the  temperature  distribution  in  both  phases. 
The  energy  equation  reads 


pC^(^^(vy)T)  =  \AT 

(2.9) 

where  C  is  the  specific  heat  at  constant  pressure,  and  X  is  the  thermal  conductivity,  both 
assumed  to  be  constant  in  the  range  of  temperatures  involved.  For  unidimensional  steady 
transport,  the  equation  (9)  reduces  to 

dT  d^T 

W—  =  K - - 

dz  dz^ 

(2.10) 

where  /c=X/pCp  is  the  fluid  thermal  diffusivity. 

Equations  (8)  and  (10)  may  be  directly  integrated,  to  give  the  distribution  of  pressure  and  of 
temperature  in  both  phases.  We  obtain 

P,  =  -P,8Z  +  q 

(2.11) 

Pi  =  -P18Z  +  C2 

(2.12) 

=  CjexpEw^z/icJ  +C4 

(2.13) 

=  Cjexpiw^z/ic,]  +Cg 

(2.14) 

as  a  function  of  the  6  integration  constants  q  to  be  determined  from  boundary  conditions. 

First  of  all,  for  z-»oo  the  temperature  in  the  vapour  phase  must  not  diverge,  which  implies 
C3=0.  Thus  the  temperature  is  constant  in  the  gas  (Tv=C4=Ti  i.e.  the  temperamre  at  the 
interface),  and  all  the  thermal  energy  is  propagated  by  the  convective  transport  of  the  latent 
heat  with  velocity  w^.  Furthermore,  we  will  neglect  the  barometric  pressure  variation  in  the 
gas  phase  (this  is  consistent  with  the  assumption  of  incompressibility  of  the  gas,  and  is 
justified  because  p,  is  small).  Thus  pv=c,=Pg,  and  the  pressure  is  constant  in  the  gas  phase. 

The  boundary  condition  at  the  rigid  conducting  plate  (z=0)  is  T=Tb,  implying  C5+C6=Tb. 
At  the  interface  z=h,  there  is  no  temperature  jump,  such  that  T,(z=h) =C5exp[W|h/fc,]  +C6=Ti. 
From  these  two  relations,  we  obtain 


^•5  = 


exp[H’/i/Kj  - 1 


^6  "  ^5 


(2.15) 
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The  constant  Cj  may  be  obtained  by  considering  the  interfacial  balance  of  normal  momentum 
(see  f.ex.  [3]),  which  reads 


Py-Pi+J'^lpv  =  0 


(2.16) 


where  the  last  term  accounts  for  the  vapour  recoil  effect,  i.e.  the  fact  that  there  is  a 
difference  of  velocity  (and  of  normal  momentum)  between  molecules  of  liquid  arriving  at  the 
interface,  and  molecules  of  vapour  leaving  it.  The  general  form  of  this  momentum  balance 
will  be  given  in  the  section  II. 3. 

The  relation  (16)  thus  gives  C2=pg+P/Pv+Pigli- 

Assuming  that  the  gas  pressure  is  given  (as  it  is  the  case  in  variable  conductance  heat  pipes), 
we  are  thus  left  with  two  unknowns  T;  (the  interfacial  temperature)  and  J  (the  mass  flux). 
Note  that  the  fluid  velocity  w,  appearing  in  eq.  15  is  linked  to  the  mass  flux  J  by  the  relation 
w,=J/p,. 

The  conservation  of  energy  at  the  interface,  the  full  form  of  which  is  given  in  the  next 
section,  expresses  that  the  jump  of  normal  heat  fluxes  is  equal  to  the  heat  used  for 
evaporation  (latent  heat  L),  plus  the  heat  transformed  into  kinetic  energy  of  the  leaving 
vapour  molecules  [3]  (note  that  the  comparatively  small  kinetic  energy  of  arriving  liquid 
molecules  is  neglected)  : 

-  X,^  .  J(L  *  *  w?)  -  /(L  *hlf)  (2.17) 

dz  dz  2  2  p^ 


This  relation  allows  to  compute  the  interfacial  temperature  T;  as  a  function  of  the  mass  flux, 
by  the  use  of  eqs  (13)  to  (15)  with  C3=0.  We  obtain 


T,  =  T,- - -exp[-vi/,WK,]) 


(2.18) 


The  last  boundary  condition,  providing  a  second  relation  between  the  mass  flux  and  the 
interfacial  temperature,  is  the  phenomenological  relation  (5).  With  the  help  of  eq.  (16),  this 
may  be  rewritten 


J 


K 


Py 


Ps^T)  -p^ 


where  1/p,  has  been  neglected  with  respect  to  l/p^ 


(2.19) 
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Relations  (18)  and  (19)  allow  us  in  principle  to  compute  T;  and  J,  provided  we  know  the 
ftmction  Ps(T),  i.e.  the  relation  between  the  equilibrium  pressure  and  the  interfacial 
temperature.  We  may  for  example  use  the  Clausius-Clapeyron  relation  [4] 

=/^oexp[-^(^-^)]  (2-20) 


where  (Po,To)  is  a  couple  of  values  lying  on  the  saturation  curve  (for  example  To=  100°C  and 
Po=l  atm  for  pure  water).  For  simplicity,  we  assume  that  the  difference  T-To  is  small,  and 
we  use  the  linearised  form  of  eq.(20) 

P,(D  =  Po*Pt(T-T„)  a-21) 


where  px=Lpo/RV  (=0.036  atm/°K  for  pure  water). 

We  may  now  proceed  to  the  resolution  of  the  system  of  equations  (18),  (19)  and  (21).  We 
first  replace  (18)  in  (21) 


Ps(T')  =  Pb  ^Pt - 

S' 


(2.22) 


where  p,,  stands  for  Po+PT(Tb-To),  i.e.  the  saturation  pressure  for  the  temperature  of  the 
bottom  plate. 

In  the  general  case,  the  analytical  resolution  of  equations  (19)  and  (22)  with  respect  to  J  and 
Tj  appears  to  be  impossible.  In  order  to  attempt  a  graphical  resolution,  we  first  solve  (19) 
with  respect  to  Ps(Ti)  : 


p  J  /2 


(2.23) 


and  both  values  of  Ps(Ti),  given  by  eqs  (22)  and  (23),  are  represented  as  a  function  of  J  in 
the  same  diagram  (fig.2).  The  solution  for  J  is  at  the  intersection  of  both  curves. 

Now,  we  may  obtain  an  analytical  result  for  J,  provided  that  its  value  is  not  too  large.  First 
of  all,  order  of  magnitude  estimates  show  that  it  is  legitimate  to  neglect  the  term  (J/pv)^  with 
respect  to  L  in  eq.  (22).  This  amounts  to  neglect  the  amount  of  kinetic  energy  imparted  to 
the  molecules  of  gas  leaving  the  interface,  compared  to  the  energy  used  for  evaporation. 
Another  simplification  is  to  neglect  the  effect  of  vapour  recoil,  i.e.  J^/p,  compared  to  other 
terms  in  (23).  This  is  valid  for  the  range  of  values  of  J  considered  in  this  analysis. 
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A  more  restrictive  approximation  concerns  the  exponential  term  in  eq. 
be  linearised  provided  that 


_  Jh  _  Wih 


<  1 


(22).  This  term  may 


(2.24) 


where  Pe  is  the  thermal  Peclet  number,  representing  the  ratio  of  the  velocity  of  the  fluid 
(associated  with  the  rate  of  phase  change)  to  the  characteristic  thermal  velocity.  It  may  also 
be  defined  as  the  ratio  of  the  thermal  relaxation  time  h^/^,  to  the  convective  time  h/w,.  Thus, 
its  value  will  be  small  provided  that  the  rate  of  phase  change  is  not  too  high.  More  precisely, 
for  a  water  layer  of  depth  h=l  mm,  this  approximation  is  valid  if  J  is  at  most  of  the  order 
10'^  g/cm^s,  representing  a  heat  transport  of  J.L=2.3  W/cm^. 


p[atm] 


Figure  2.2:  Graphical  determination  of  J.  Pg=l  atm  is  the  vapour  pressure,  Pb=1.5  atm  is 
the  saturation  pressure  at  the  rigid  plate  temperature  T|,=487K  (eq.  21).  Curves  1  and  2  are 
given  by  equations  (22)  and  (23),  curve  3  is  the  Clausius-Clapeyron  curve  (eq.20)  in  which 
the  temperature  of  the  interface  is  evaluated  by  eq.  18.  All  the  parameters  are  those  of  water 
(T(,=473K,  Po=1  atm),  with  a  depth  h=lmm,  and  an  accomodation  coefficient  /3= 0.1.  The 
supersaturation  is  Tb-To=14°.  It  is  found  that  1=2.6  10"^  g/cm^s  (J.L=0.6  W/cm^). 


With  these  approximations,  eq.  (22)  may  be  written 


■'pi 


p.K 


(2.25) 
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and  eq.  (23)  reads 

pm  ‘P,*^ 

Equalisation  of  these  two  relations  leads  to  the  mass  flux 

.  _  P»  -Ps 

P,  PiLh 

This  relation  involves  only  measurable  quantities,  at  the  exception  of  the  phenomenological 
coefficient  K.  However,  as  mentioned  earlier,  K  can  be  roughly  estimated  by  the  formula 
K=Pp,(T)  (M/2'kRT)^'^.  For  water,  this  gives  K=6  lO'*®  gs/cm'*  (with  a  value  0.1  of  the 
accommodation  coefficient  |S).  Then,  sincepxLh/X|=:8  10*cm/s>  >Pv/K=10®cm/s,  the  first 
term  in  the  denominator  of  (27)  is  negligible  compared  to  the  second  one  (quasi-equilibrium 
case),  and  we  get  the  relation 


(2.26) 


(2.27) 


j  =  (2.28) 

PjLh 

However,  in  the  following,  we  will  adopt  the  expression  (27)  of  the  reference  mass  flux  J, 
in  part  because  the  actual  values  of  K  could  be  quite  different  from  those  estimated  with  the 
help  of  the  Hertz-Knudsen  relation.  Moreover,  the  value  of  the  accommodation  coefficient 
P  is  generally  found  to  be  significantly  reduced  by  the  presence  of  impurities  on  the 
interface,  such  that  it  could  not  be  legitimate  to  neglect  non-equilibrium  effects  (finite  K). 

We  now  mm  to  a  stability  analysis  of  the  reference  state  investigated  in  this  section. 


II.3.  Stability  of  the  reference  solution  versus  hydrodynamic  fluctuations 

We  now  have  to  remm  to  a  more  general  formulation  of  the  problem,  which  in  particular 
allows  for  two-dimensional  flucmations  of  velocities,  of  pressure  and  of  ternperamre,  and 
involves  surface  tension  effects.  As  explained  in  the  introduction,  surface  tension  gradients 
(caused  by  gradients  of  ternperamre)  might  be  able  to  destabilise  the  above-calculated 
solution,  provided  that  the  ternperamre  gradient  perpendicular  to  the  interface  (created  by 
evaporation)  exceeds  a  critical  value. 

In  this  section,  we  define  the  linear  stability  problem  to  be  smdied.  The  system  of  equations 
and  boundary  conditions  will  be  written  under  a  dimensionless  form,  and  a  particular  model 
of  the  problem  will  be  introduced.  This  is  known  in  the  literamre  as  the  one-sided  model  [3], 
for  which  the  dynamics  of  the  gas  phase  can  be  decoupled  from  the  dynamics  in  the  liquid 
phase.  We  may  thus  solve  the  problem  (with  suitable  boundary  conditions)  in  the  liquid  phase 
only.  This  model  requires  the  assumption  that  the  ratios  of  volumic  masses  p^/p,,  of  thermal 
conductivities  \/\  and  of  dynamic  viscosities  are  small  compared  to  unity,  which  is 
verified  for  most  liquids  in  contact  with  their  vapour,  far  enough  from  the  critical  point.  It 
will  be  demonstrated  that  the  effect  of  evaporation  can  be  modeled  (in  a  first  approximation) 
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by  introducing  a  free  surface  dimensionless  Biot  number  Bi  (for  which  a  value  is  derived  in 
this  section). 

In  section  II.4,  the  linear  stability  problem  is  solved,  and  thresholds  for  the  Marangoni- 
Benard  instability  are  provided.  The  obtained  solutions  are  valid  as  long  as  the  amplitudes 
of  perturbations  around  the  reference  solution  are  small.  In  section  III,  the  full  nonlinear 
Marangoni-Benard  problem  will  be  investigated  (finite-amplitude  solutions)  for  short  waves 
(perturbations  located  near  the  interface),  after  having  simplified  the  problem  by  considering 
the  case  of  an  infinite  Prandtl  number  fluid  (which  leads  to  results  applicable  to  the  case  of 
normal  liquids,  including  water  for  which  Pr=7). 

i)  Model  equations 

The  method  used  to  study  the  stability  of  the  reference  solution  (computed  in  the  previous 
section)  is  general,  and  proceeds  from  the  definition  of  stability,  i.e.  the  resistance  of  the 
basic  solution  to  infinitesimal  hydrodynamic  fluctuations.  The  equations  governing  the 

evolution  of  infinitesimal  perturbations  V ,  p’  and  T’  of  velocity,  pressure  ai^  temperatore 
around  their  reference  state  values  and  are  obtained  by  inserting  v  =  +  v' , 

p=p,^+p’  and  T=Tref+T’  in  the  incompressibility  relation  7.7  =  0  and  in  eqs  (7)  and  (9), 
and  by  remarking  that  the  reference  solution  satisfies  these  equations.  We  then  obtain 


V.v  =  0 

(2.29) 

+pw^^w  =  pAw  -Dp 

(2.30) 

p  —  +  pwDu  =  pAu  -  ^ 
dt  dx 

(2.31) 

kAT 

(2.32) 

where  the  primes  have  been  omitted  for  simplicity,  all  second-order  terms  involving  products 
of  perturbations  have  been  neglected,  both  components  of  the  Navier-Stokes  equation  (7) 

have  been  separated  (7=  +  wT^),  D  denotes  the  z-derivative,  and  K=X/pCp  is  the  thermal 

diffusivity.  Note  that  in  the  following,  subscripts  v  and  1  will  be  omitted  when  referring  to 
a  liquid  quantity. 

On  the  rigid  conducting  bottom  plate  z=0.  The  boundary  conditions  for  the  perturbations  are 

u  =  w  =  T  =  0 

because  the  total  temperature  is  maintained  constant  at  Tb  (no  perturbation)  and  the  total 
velocity  is  kept  at  its  reference  value  V  =  (Jref  is  the  reference  value  of  the 
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mass  flux,  given  by  eq.  27). 

At  the  interface  z=h+^(x,t),  the  linearisation  of  the  mass  conservation  relation  (4)  leads  to 

As  Jref==PiWref,  the  peiturhed  mass  flux  J’  is  given  by 

J'  =  p^(w,(z=h)-w^)  (2-33a) 

and  will  be  denoted  by  J  in  the  following.  Note  that  an  equivalent  treatment  of  the  second 
equality  in  (4)  gives 


J'  =  Py(wXz=h)-w^) 


(2.33b) 


The  general  form  of  the  momentum  balance  at  the  interface  reads  [1] 


Ip((v-Vj;).«)  v+pn  +  T.n} 


odtldx+  tda!  dx 
_ 


(2.34) 


where  I  x  I  stands  for  the  discontinuity  Xv-x,  (v  and  1  denote  respectively  the  vapour  and  the 
liquid),  '«=(-^M^  +T^)/iVand  T=(T^  +  are  respectively  the  free  surface  normal  and 

tangential  unit  vectors  (N=(l+^’^)‘^^  is  the  normalisation  factor),  a  is  the  surface  tension,  t 
is  the  viscous  stress  tensor  with  components  Ti^=-p(dVi/dx^+dv/dx). 


Applying  the  one-sided  model  approximation  (pJpiM),  pJpfM)),  substituting  p^pref+P. 
w-^w,ef+w,  and  linearising,  the  projection  of  (34)  on  the  surface  normal  reads 


-  ^  -  J, 


(PXef  -  (P)r^+Py  'Pi  =  (t,  .  «). /I  -  -  2 +  U 


Jr,J 


32^ 

3x2 


(2.35) 


where  the  relation  (33)  has  been  used.  This  relation  must  be  expressed  at  z=h-l-|(x,t).  Since 
^(x,t)  <  <  1,  the  Taylor  series  around  z=h  can  be  limited  to  linear  terms  in  ^  Using  (16) 
at  z=h,  we  obtain 


Py-Pt 


-p^g^-2ppw,-2^ 


+  a 


3x2 


(2.36) 


We  must  also  project  the  relation  (34)  on  the  surface  unit  tangent  t.  This  gives,  after 
linearisation. 


du,  3w,  ^ 
u,(—+ — -)  =  O. 
dzdx 


^  3x 


(2.37) 


where  aj  is  the  surface  tension  variation  with  temperature.  We  may  use  the  continuity 
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relation  (29)  to  rewrite  (37)  under  the  form  (expressed  at  z=h) 


,  d^w  s  d'^T 

■  a?’  ■  aF 


Note  that  in  deriving  eq.  (37),  we  used  the  no-slip  condition 

iT  =  v^.7 

whose  perturbed  and  linearised  form  is 


(2.39) 


=  (2.40) 

dx 

Now,  perturbing  eq.  15,  and  applying  the  one-sided  model  approximations,  we  get  the 
perturbed  mass  flux  as 


J  =  K 


Pv  Pi 


El 

Pv 


(2.41) 


We  may  eliminate  p^  between  eqs  (36)  and  (39),  and  get  the  following  relation  for  J  ; 


J  = 


K 

2KJ, 

T 

Pv 


Ei(t+^dt  )-Ei  + 

f>V  Pv 


Ell^  +  2  ELDw  I 


Pv 


(2.42) 


Having  expressed  the  conservation  of  mass  and  of  momentum  at  the  interface,  we  still  have 
to  write  the  general  conservation  equation  for  the  energy  at  the  interface,  which  reads  (in  the 
one-sided  model)  as  [3] 


JiL  +  h-?)  =  -\n.VT, 

2  Pv 


(2.43) 


The  perturbed  and  linearised  form  is  derived  in  the  usual  way.  We  obtain 

-\DT,  -  =  J(L  +1(:^)2) 

"  Pv 


(2.44) 


ii)  Dimensionless  equations 

The  equations  and  boundary  conditions  can  be  put  under  dimensionless  form  by  using  h  (the 
unperturbed  liquid  depth)  as  unit  length,  h^/zc  as  unit  time,  zc/h  as  unit  velocity,  and  /xzc/h^  as 
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unit  pressure  (fluid  quantities  without  indices  are  those  of  the  liquid  phase).  The  temperature 
unit  d  is  based  on  the  amplitude  of  the  thermal  gradient  at  the  interface  in  the  basic  state 
(found  from  eq.  14,15,18  and  24)  : 


e  -h\  DT^(z=h)  I - - — 

The  dimensionless  equations  valid  in  the  liquid  phase  thus  read  (the  same  notations  are  kept 
for  dimensionless  velocity,  temperature  and  pressure)  : 


Pr(Aw  -  Dp)  ■■ 

Pr(Au-^) 

ox 

Ar  +  wexp[Pe(z-] 


0 

(2.45) 

—  +PeDw 

(2.46) 

dt 

—  +  PeDu 

(2.47) 

dt 

\  =  ^  +  PeDT 

(2.49) 

where  Px^filpK  is  the  Prandtl  number  of  the  liquid,  and  the  Peclet  number  Pe=w,efh//c  has 
been  defined  earlier  (eq.  24). 

The  boundary  conditions  are  : 

-  At  the  bottom  plate  (z=0)  : 


u  =  w  =  T  =  0 


(2.50) 


-  At  the  interface  tz=l)  ; 

The  mass  conservation  relation  -  definition  of  mass  flux  (eqs.  33a  and  42)  : 


J  =  Pr'*(iv-Wj.)  =  Pr'^PeBi 


(T-^)-U,(p-2Dw)+Il^^-n, 


(2.51) 


In  this  relation,  the  Biot  number  Bi  has  been  defined  by 

Bi  =  ^ 

1-2K^ 


(2.52) 


where  (J^f/pv)^  has  been  neglected  with  respect  to  L  in  the  expression  of  d  (as  for  eqs  25-28). 
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Other  dimensionless  groups  appearing  in  eq.  51  are  : 


l^h^PrL 


Pi8\ 

'^rifPr^ 


Hb  = 


a\ 


J^^h^PrL 


Another  boundary  condition  at  z=l  is  the  normal  momentum  balance  (eq.  36) 

p  -p,+2Dw,  =  -Ga^ 

'  dx^  p*Pr 

where  we  introduced  the  density  ratio  p*=pjp\,  the  Crispation  number 


(2.53) 


which  can  be  viewed  as  a  measure  of  the  free  surface  deformability,  and  the  Galileo  number 

Vi 


characterising  the 
We  also  have  the 


influence  of  gravity  on  surface  deformations, 
tangential  stress  condition  (eq.  38) 


dx^  dx^ 


(2.54) 


where  the  Marangoni  number,  characterising  the  destabilising  influence  of  the 
thermocapillary  effect,  is  defined  by 

-Oj-dh  -Oj.^h^ 

Ma  =  -J—  =  — - 

pK  pK 

where  13  =  \  DT,^f(z=h)  !  is  the  temperature  gradient  amplitude  at  the  interface. 
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The  no-slip  condition  (eq.  40)  reads 


(2.55) 


Ui  =  u^  + 


D*  dx 


The  thermal  boundary  condition  is  obtained  from  eq.  44,  by  remarking  that  D2T„f=-Pe0/h2. 
Using  (42)  and  the  definitions  of  Bi  and  11;,  i= 1,2,3,  we  get 


DT-Pe^  =  Bi 


+Il,(p-2Dw)-U^^  +U, 


32^ 

3x2 


(2.56) 


and  the  last  boundary  condition  is  the  kinematic  condition,  whose  dimensionless  linearised 
form  reads 


dt 


(2.57) 


The  system  of  equations  and  boundary  conditions  (45)-(57)  is  still  not  closed,  since  some 
quantities  referring  to  the  vapour  phase  still  appear  in  boundary  conditions  (53)  and  (55). 
However,  an  important  limiting  case  occurs  when  the  Crispation  number  is  very  small  (for 
a  water  layer  of  depth  h=lmm,  Cr=1.4  10®),  such  that  ^->0  according  to  the  boundary 
condition  (53).  It  will  be  shown  in  the  next  section  that  this  approximation  is  valid  for 
perturbations  with  sufficiently  small  wavelengths,  for  which  surface  tension  strongly  damps 
surface  deformations.  Note  that  for  large  wavelengths,  gravity  is  expected  to  play  an 
equivalent  stabilising  role  (the  Galileo  number  Ga=  10®  for  1mm  water).  Order  of  magnitude 
estimates  also  show  that  Hi  =  lO'^-lO-^  for  the  range  of  mass  fluxes  considered,  such  that  the 
pressure  term  may  be  neglected  in  eqs.  51  and  56,  provided  |  p-2Dw  |  / 1 T  |  is  not  larger  than 
10'®  (this  will  be  verified  a  posteriori). 

With  these  approximations,  the  system  of  equations  (45)  to  (49)  is  unchanged,  while  the  free 
surface  boundary  conditions  51-57  reduce  to  : 


w  =  PeBi  T 

(2.58) 

D^w  - -  =  Ma - 

(2.59) 

3x2 

DT  +  BiT  =  0 

(2.60) 

and  the  system  is  closed  in  this  limit,  since  no  quantities  linked  to  the  gas  phase  appear  in 
these  boundary  conditions. 

It  was  also  assumed  in  the  last  section  (eq.  24)  that  the  liquid  Peclet  number  Pe  is  small 
compared  to  unity.  In  the  limit  Pe^,  the  full  problem  just  derived  reduces  to  the  problem 
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linearly  treated  by  Pearson  [5]  (see  next  section).  However,  in  Pearson’s  analysis,  the  value 
of  the  Biot  number  Bi  is  not  related  to  evaporation  parameters  (because  the  heat  transfer  at 
the  free  surface  is  described  by  a  classical  Newton’s  cooling  law,  with  a  phenomenological 
coefficient  whose  estimation  may  be  expected  to  depend  on  experimental  peculiarities,  such 
as  the  gas  depth,  thermal  conductivity,  gas  convection,  ...).  In  our  analysis,  the  relation  (52) 
will  be  used  to  evaluate  the  Biot  number.  In  spite  of  the  fact  that  this  relation  also  contains 
a  phenomenological  coefficient  K  (defined  by  eq.  5),  it  can  be  expected  that  a  satisfactory 
order  of  magnitude  for  K  can  be  obtained  from  the  kinetic  theory  of  gases  (for  water  with 
h=lmm,  and  an  accommodation  coefficient  /8=0. 1,  we  find  Bi=:  10^,  but  this  value  could 
be  significantly  reduced  by  impurities  leading  to  much  lower  values  of  the  accommodation 
coefficient  /3=0.01  or  even  ,8=0.001),  and  that  this  value  will  be  less  dependent  on 
convection  in  the  gas  phase  than  it  is  for  the  Newton’s  heat  transfer  coefficient.  Furthermore, 
the  Biot  number  computed  from  eq.  52  will  prove  to  be  useful  for  investigating  the  influence 
of  non-equilibrium  effects  (finite  K). 

II.4.  Normal  modes  analysis 

In  the  limit  Pe-^,  which  can  be  eonsidered  as  a  good  approximation  for  reasonable  mass 
fluxes  and  small  liquid  depths,  the  system  of  eqs  (45)  to  (49)  reduces  to 


V.v  =  0 

(2.61) 

> 

1 

ii 

1 

(2.62) 

(2.63) 

dx  dt 

A-r 

AT  +  w  =  — 

(2.64) 

dt 

with  the  boundary  conditions  u=w=T=0  at  z=0.  At  the  interface  z=h 

[,  from  eqs  58-60, 

w  =  0 

(2.65) 

a2T 

D^w  =  Ma-— 

(2.66) 

dx^ 

DT  +  BiT  =  0 

(2.67) 

This  problem  has  been  solved  by  Pearson  [5]  in  the  neutral  stability  case  (8/3t=0),  and  the 
general  resolution  will  now  be  presented. 

It  is  possible  to  reduce  the  problem  (61)-(67)  to  a  problem  for  w  and  T  only.  This  is  done 
by  deriving  eq.  63  with  respect  to  x,  and  adding  the  result  to  the  z-derivative  of  eq.  62. 
Using  eq.  61,  we  obtain  Ap=0.  Now,  applying  the  Laplacian  operator  A  to  eq.  62,  we 
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obtain  a  fourth  order  equation  for  w  : 


Aw 

dt 


(2.68) 


Thus,  we  only  need  to  solve  eqs  68  with  eq.  64,  and  apply  the  boundary  conditions 
w=Dw=T=0  at  z=0  (the  condition  Dw=0  comes  from  u=0  at  z=0  and  eq.  61),  and 
interfacial  conditions  65-67. 

Separation  of  variables  shows  that  the  general  solution  of  the  linear  problem  is  the 
superposition  of  Fourier  modes  that  may  be  written 


=  exp[(r^r  +  ikx] 


w,(z) 

w 


(2.69) 


where  the  parameter  k  is  the  wavenumber  of  the  disturbance  (the  wavelength  is  27r/k). 


By  introducing  eq.  69  in  the  above  problem,  it  is  seen  that  each  Fourier  component  (normal 
mode)  must  satisfy  the  eigenvalue  problem 


(2.70) 


(p2-k^)T,(z)  +w^iz)  =  (Tji^iz) 


(2.71) 


with  boundary  conditions 

w,(0)  =  Dw,(0)  =  r,(0)  =  0  (2-72) 

w/l)  =  +  k^MaT,(l)  =  DT,(l)  +BiT,il)  =  0  (2-73) 


It  is  possible  (although  calculations  are  not  reproduced  here  for  simplicity)  to  write  the 
general  solution  of  the  ordinary  differential  equations  70-71,  which  depends  on  6  constants 
to  be  determined  from  the  boundary  conditions  72-73.  This  provides  a  system  of  6 
homogeneous  equations  for  the  6  unknown  coefficients,  which  admits  a  non-trivial  solution 
if  and  only  if  a  compatibility  condition  (cancelation  of  the  determinant  of  the  matrix  of  the 
system)  is  satisfied.  This  characteristic  relation  may  be  written 

Aia^,k,Pr,Bi,Ma)  =  0  (2-74) 

and  allows  to  compute  the  growth  constants  for  given  k,  Bi  and  Ma.  Then,  the  reference 
state  will  be  stable  provided  all  growth  constants  have  negative  real  part  for  all  k.  At  the 
contrary,  it  will  be  unstable  if  there  exist  some  values  of  k  for  which  at  least  one  eigenvalue 
has  a  positive  real  part  (from  eq.  69,  it  is  seen  that  it  will  grow  exponentially  in  time).  In 
practice,  we  determine  the  threshold  of  the  instability  by  substituting  (the  possibility 
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of  oscillatory  onset  a=iw  can  be  rejected,  as  shown  in  [6])  in  eq.  74.  After  simplifying  the 
so-called  neutral  stability  relation,  we  obtain 

_  \6k  {kcoshjk)  +Bi sinh(^))  (sinh(2 A:)  - 2k)  (2.75) 

sinh(3A:)  -  3  sinh(^)  -Ak'^  cosh(^) 

which  is  the  Pearson’s  result  [5].  Note  that  this  result  does  not  depend  on  the  Prandtl  number 
Pr  (when  (rk=0.  the  problem  70-73  does  not  depend  on  Pr  anymore).  This  relation  is 
represented  on  fig.  3  for  several  values  of  the  Biot  number  Bi.  For  given  Ma  and  Bi,  it  can 
be  verified  that  the  reference  state  (computed  in  section  II. 2)  is  stable  against  a  disturbance 
with  wavenumber  k,  provided  the  point  (k,Ma)  lies  below  the  curve  for  the  given  Biot 
number  Bi. 


Ma 


Figure  2.3  :  Neutral  stability  curves  for  various  Biot  numbers  (indicated  for  each  curve),  as 
given  by  eq.  75.  The  critical  Marangoni  number  is  given  by  the  minimum  these  curves  as 
a  function  of  the  wavenumber  k  (dashed  line). 


The  general  solution  of  the  stability  problem  may  be  written  as  the  superposition 


w(.x,z,f) 

T(x,z,t) 


I  expfff^t  +  ikx] 


'w,(z) 


(2.76) 


with  arbitrary  Fourier  coefficients  a^.  Then,  the  critical  Marangoni  number  Ma,  above  which 
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convection  sets  in  is  detennined  by  the  minimum  of  the  neutral  stability  curve  eq.  75.  For 
Bi=0  for  example,  it  is  found  that  Mae=79.6  for  k,=1.993.  This  value  is  increasing  with 
the  Biot  number,  as  seen  on  figs  3  and  4. 

Just  above  the  critical  point  Ma=Mac,  disturbances  with  wavenumbers  lying  in  a  very  narrow 
range  centred  around  k=kc  will  be  amplified,  while  disturbances  with  wavenumbers  outside 
this  range  will  be  damped.  Thus,  after  some  time,  the  liquid  layer  is  in  a  convection  regime, 
with  a  structure  that  looks  like  the  critical  mode  with  a  unique  wavenumber  k^.  Figure  5 
represents  this  convection  structure,  obtained  by  multiplying  the  critical  mode 


TAz,t) 


=  exp[/A:jc] 


w,(z) 

TAz) 


+  c.c. 


by  a  reasonable  amplitude  a^e-  In  reality,  the  choice  of  this  amplitude  is  for  the  moment 
arbitrary,  because  the  linear  analysis  just  performed  does  not  allow  its  calculation. 


Note  that  the  basic  temperature  profile  (-z  in  dimensionless  form)  has  been  added  to  the 
temperature  disturbance  in  representing  the  isotherms  of  fig.  5. 


II.5.  Importance  of  a  nonlinear  analysis 

The  linear  analysis  just  performed  does  not  allow  to  determine  the  amplitudes  a|j  in  the 
Fourier  decomposition  76.  In  fact,  this  is  not  the  only  limitation  of  linear  analyses.  For 
example,  it  is  seen  from  eq.  76  that  once  the  threshold  of  instability  is  exceeded,  some  modes 
with  positive  growth  rate  will  grow  exponentially  in  time,  such  that  their  amplimde  can  be 
expected  to  grow  indefinitely.  This  unrealistic  result  is  due  to  the  fact  that  the  perturbations 
have  been  assumed  infinitesimal,  such  that  nonlinearities  of  the  basic  equations  have  been 
neglected  (see  section  II.3,  eqs  30-32).  When  Ma>Mae,  some  perturbations  grow  and  may 
no  longer  be  assumed  infinitesimal,  such  that  convective  nonlinearities  have  to  be 
reincorporated  in  the  balance  equations  (see  eqs  7  and  9).  This  will  be  done  in  the  next 
section,  and  it  will  be  seen  that  the  amplitudes  3^(0  will  saturate  to  some  finite-amplimde 
values. 

These  nonlinear  terms  will  also  lead  to  couplings  between  the  Fourier  modes  with  different 
wavenumbers  k,  which  are  decoupled  at  the  linear  stage.  More  generally,  the  linear  problem 
might  have  been  solved  in  three  dimensions  (one  vertical  direction  z,  and  two  horizontal 
directions  x  and  y),  and  we  would  have  obtained  the  same  equations  64  and  68,  with 
boundary  conditions  65  to  67  where  d^/dx^  has  to  be  replaced  by  the  horizontal  Laplacian 
operator  A,=d^/dx'^+d^/dy^.  Then,  the  Fourier  integral  76  would  have  been  replaced  by 


w(x,z,t) 

T(x,z,t) 


dk^dky  a^QX^[aA  +  i(k^x+k^y)] 


wAz) 

m 


(2.77) 


However,  due  to  the  isotropy  in  the  horizontal  directions,  we  would  have  obtained  the  same 
spectral  problem  70-73,  with  the  wavenumber  k=(k,2+ky2)i/2.  Then,  the  linear  problem  leads 
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Mae 

1500 


Figure  2.4  :  Critical  Marangoni  number  Ma<.  (top)  and  critical  wavenumber  (bottom)  as 
a  function  of  the  Biot  number  Bi.  For  Bi->oo,  Mac=32.1  Bi  and  kc=3.01. 
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Figure  2.5:  Critical  perturbations  (above  :  streamlines,  below:  isotherms).  Bi=0  (kc=2, 
leading  to  periodic  convection  with  wavelength  X=27r/kc=:3.14  in  units  of  the  layer  depth). 


to  the  same  characteristic  relation,  and  at  the  critical  point,  not  only  one  but  all  the  modes 
lying  on  the  critical  circle  gj-g  critical  (the  problem  is  thus  infinitely 

degenerated).  The  evolution  of  the  amplitudes  a(k^,ky)  of  these  modes  is  decoupled  at  the 
linear  stage  (all  grow  in  time  with  the  same  growth  constant),  but  as  soon  as  they  leave  the 
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Figure  2.6  :  Top  view  of  the  free  surface  isotherms,  for  3  roll  patterns  (a,b,c)  with 
orientations  differing  by  60°.  Their  superposition  with  equal  amplitudes  leads  to  hexagons 
(d).  Black  regions  :  cold,  white  regions  :  hot. 
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linear  domain,  nonlinear  interactions  between  them  will  provide  couplings  (and  competition) 
of  Fourier  modes  (convection  rolls  with  axes  perpendicular  to  the  wavevector 

T  =  with  different  orientations.  In  section  III,  the  competition  between  3  roll 

patterns  with  wavevectors  lying  on  the  critical  circle  and  forming  angles  of  60°  with  each 
other  will  be  investigated,  as  it  is  known  that  this  leads  to  the  experimentally  observed 
hexagonal  structures  (see  fig.  6). 

In  addition  to  this  orientational  degeneracy  of  the  linear  problem,  there  exist  another 
"bandwidth"  degeneracy,  linked  to  the  fact  that  at  a  given  supercritical  Marangoni  number 
Ma,  a  band  of  disturbances  is  amplified  (all  the  wavenumbers  lying  above  the  neutral  stability 
curve).  In  this  band,  there  exist  a  wavenumber  possessing  the  maximal  growth  rate.  Thus, 
at  the  linear  stage,  this  mode  will  be  amplified  more  quickly,  and  will  dominate  other 
perturbations.  Here  again,  nonlinear  couplings  between  modes  with  different  k  will  occur, 
and  some  disturbances  will  develop  at  the  expense  of  others,  and  determine  the  wavelength 
of  the  final  convective  structure,  which  may  then  not  be  determined  from  the  linear  analysis. 
This  wavelength  selection  problem  will  now  be  examined  for  two-dimensional  perturbations, 
and  assuming  that  the  layer  of  liquid  is  sufficiently  deep  (short- wavelength  perturbations). 
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II.  6.  Conclusions  of  section  II 


In  this  section,  we  developed  a  model  describing  evaporation  of  a  liquid  layer  in  contact  with 
its  vapour,  and  lying  on  a  heated  plate  maintained  at  constant  temperature.  For  a  given 
pressure  Pg  of  the  gas  phase,  evaporation  occurs  if  the  temperature  of  the  bottom  plate  is 
superior  to  the  saturation  temperature  corresponding  to  the  pressure  Pg  (obtained  from  the 
classical  Clausius-Clapeyron  equation).  If  the  temperature  of  the  bottom  plate  is  inferior  to 
this  temperature,  then  condensation  occurs.  In  order  for  the  system  to  reach  a  steady 
(reference)  state,  we  assume  that  fresh  liquid  is  injected  at  the  bottom  plate,  at  a  rate  which 
exactly  matches  the  rate  of  evaporation  (or  condensation)  at  the  free  surface.  This  is  done 
in  order  to  mimic  steady  states  reached  in  evaporating  devices  such  as  heat  pipes.  A 
graphical  determination  of  the  steady  mass  flux  is  proposed,  together  with  an  analytical  result 
based  on  a  linear  approximation  valid  for  mass  fluxes  that  are  not  too  large  (more  rigorously, 
the  Peclet  number  of  the  liquid  needs  to  be  small  compared  to  unity). 

Our  model  of  evaporation  includes  the  possibility  of  non-equilibrium  effects,  i.e.  the  interface 
is  not  assumed  to  be  in  a  thermodynamic  equilibrium  state,  for  which  chemical  potentials  on 
each  side  of  the  interface  are  equal  (the  temperature  of  the  interface  would  then  be  equal  to 
the  saturation  temperature  computed  from  the  Clausius-Clapeyron  equation).  Instead,  the 
difference  of  these  chemical  potentials  is  considered  as  a  generalised  force  (in  the  sense  of 
irreversible  thermodynamics)  that  generates  a  mass  flux  J  proportional  to  it  (in  the  linear 
approximation).  The  proportionality  coefficient  K  is  a  phenomenological  coefficient  similar 
to  a  diffusion  coefficient,  or  a  thermal  conductivity,  and  its  value  is  not  presently  well 
known.  By  comparing  the  expression  we  obtain  for  J  with  the  classical  Hertz-Knudsen 
equation,  rough  estimations  are  obtained  for  K,  which  however  depend  on  the  value  of  the 
accommodation  coefficient  (3,  for  which  large  discrepancies  exist  between  theory  and 
experiment.  We  thus  considered  the  coefficient  K  as  a  free  parameter,  the  limit  K-^oo 
defining  the  state  of  thermodynamic  equilibrium  of  the  interface,  and  the  limit  K=0 
describing  the  case  of  an  "impermeable"  interface  (J=0). 

In  a  second  stage,  the  analysis  of  the  hydrodynamic  stability  of  the  computed  reference  state 
was  performed,  taking  into  account  surface  tension  variation  with  temperature.  A  system  of 
equations  and  boundary  conditions  describing  a  one-sided  model  of  convection  (the  dynamics 
of  the  liquid  layer  is  decoupled  from  the  dynamics  of  the  gas)  was  developed  and  put  under 
dimensionless  form.  Order  of  magnitudes  estimates  of  the  dimensionless  numbers  entering 
the  problem  shows  that  some  effects  can  be  neglected.  These  are  the  surface  deformation  (the 
crispation  number  Cr  is  vanishingly  small),  and  the  nonlinearity  of  the  basic  temperature 
profile  created  by  evaporation  (the  Peclet  number  is  small  compared  to  unity,  especially  for 
small  liquid  depth).  It  is  finally  shown  that  the  problem  of  stability  reduces  to  the  case  of  a 
motionless  liquid  layer  in  contact  with  an  inert  gas  phase  at  an  undeformable  interface  whose 
surface  tension  varies  with  temperature,  the  basic  temperature  profile  in  the  layer  being  linear 
(Pearson's  problem).  Evaporation  is  incorporated  via  the  presence  of  a  free  surface  heat 
transfer  coefficient  (a  Biot  number),  for  which  we  have  obtained  a  useful  formula  depending 
on  the  characteristics  of  the  liquid  (latent  heat,  Clausius-Clapeyron  slope,  thermal 
conductivity  and  density)  and  on  the  phenomenological  coefficient  K.  As  the  critical 
Marangoni  number  (defining  the  critical  value  of  the  thermal  gradient,  or  equivalently  of  the 
evaporation  mass  flux)  increases  with  the  Biot  number  (which  increases  with  K),  it  can  be 
concluded  that  evaporation  stabilises  the  layer  with  respect  to  Marangoni  convection. 
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III.  FINITE  AMPLITUDE  REGIMES  OF  THE  SHORT-WAVE  INSTABILITY 


In  this  section,  a  model  of  the  pure  thermocapillary  instability  in  layers  of  infinite  depth  is 
developed  in  the  framework  of  the  amplimde  equations  formalism.  We  make  use  of  the 
eigenfunctions  at  a  given  Marangoni  number  Ma  (as  determined  in  the  previous  section)  as 
a  basis  for  the  nonlinear  problem,  rather  than  the  neutral  stability  functions  (as  it  is  often 
done  for  weakly  nonlinear  analyses).  It  will  be  shown  that  third  order  equations  may  visibly 
be  extrapolated  rather  far  above  the  threshold.  In  particular,  results  will  be  obtained  about 
the  wavelength  selection  problem  between  fastest  growing  modes  (wavenumbers  around 
for  a  zero  free  surface  Biot  number)  and  critical  modes  {k^-^O  and  Ma^-^0  when 
the  layer  depth  is  infinite).  Transient  numerical  integration  of  the  equations  reveals  an 
unbounded  growth  of  the  mean  wavelength,  thus  indicating  the  absence  of  an  intrinsic 
wavelength  for  this  physical  system.  This  is  explained  in  terms  of  the  mean  (horizontally 
averaged)  temperature  profile  distortion  by  convection.  The  final  steady  state  of  this  evolution 
(imposed  wavelength)  is  then  approximated  analytically.  Earlier  results  about  the  competition 
between  rolls  and  hexagonal  patterns  are  qualitatively  recovered.  These  solutions  are  then 
investigated  in  the  limit  Ma-^oo,  where  power  law  relationships  are  derived  for  main 
convective  quantities.  In  particular,  a  saturation  behaviour  is  obtained  for  a  quantity  (the  bulk 
temperature  decrease),  which  can  be  considered  as  a  measure  of  the  heat  transport  increase 
due  to  convection. 

The  contents  of  this  section  form  the  subject  of  one  publication  entitled  "Finite-amplitude 
regimes  of  the  short-wave  Marangoni-Benard  convective  instability",  P.  Colinet,  J.C.  Legros, 
Y.  Kamotani,  P.C.  Dauby,  and  G.  Lebon,  to  appear  in  Phys.  Rev.  E,  August  95. 

III.l.  Introduction 

When  the  threshold  of  the  Marangoni-Benard  instability  is  exceeded,  various  dissipative 
structures  are  experimentally  observed,  some  of  them  localised  near  the  interface,  in  the  form 
of  small  cells  eventually  embedded  in  larger  convective  structures  [5],  soliton-like 
propagating  waves  [6],  or  interfacial  turbulence  [7].  In  other  conditions,  Marangoni-Benard 
convection  is  rather  similar  to  the  classical  buoyancy  induced  Rayleigh-Benard  convection 
[8-11],  with  patterns  extending  far  into  the  bulk  of  surrounding  liquids,  and  actually  reaching 
(and  influenced  by)  the  boundaries  of  the  experimental  vessel.  Note  that  these  apparently 
different  forms  of  convection  generally  result  in  substantial  increases  of  heat/mass  transfer 
through  the  interface. 

Figure  1  reproduces  the  neutral  stability  results  of  Pearson  [13],  as  well  as  results  obtained 
by  Scanlon  and  Segel  [14]  in  the  case  of  a  layer  of  infinite  depth  (the  length  d  appearing  in 
Ma=-o-p(3d^/fXK  then  represents  an  arbitrary  length).  At  a  given  Ma,  perturbations  with 
wavenumbers  in  the  range  lying  above  the  neutral  stability  curve  possess  a  positive 
amplification  rate.  There  also  exists  a  particular  wavenumber  in  this  range  which 
possesses  the  maximal  amplification  rate.  In  slightly  supercritical  conditions  {Ma=^Ma^,  k„,^ 
is  close  to  the  critical  wavenumber  k^,  and  does  indeed  predict  the  size  of  convective 
(hexagonal)  cells  observed  experimentally.  When  Ma  is  increased,  k„,^  is  seen  to  increase, 
and  the  prediction  of  the  selected  wavelength  becomes  quite  complicated,  since  it  involves 
non-linear  competition  between  modes  in  the  unstable  range.  Progress  has  recently  been 
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is  not  the  case  for  the  Marangoni-Benard  instability  (the  critical  wavenumber  in  a  finite  layer 
of  depth  h  scales  as  1/h,  such  as  for  the  Rayleigh-Benard  instability).  However,  we  may 
conjecture  that  this  does  not  rule  out  the  possibility  of  an  intrinsic  wavelength,  linked  to  the 
presence  of  fastest  growing  modes  (as  it  seems  to  be  the  case  for  the  Rayleigh-Taylor 
instability  [19]).  Note  that  the  finite  wavenumber  of  the  fastest  growing  mode  generally 
depends  on  the  driving  force  amplitude  (the  thermal  gradient  in  our  case).  Attempts  to 
answer  to  the  question  of  the  preference  of  such  modes  at  a  given  supercritical  driving  force 
obviously  have  to  incorporate  nonlinear  effects  in  the  analysis. 

After  having  described  the  model  in  section  III. 2,  in  which  the  derivation  of  weakly  nonlinear 
results  is  also  described,  we  discuss  this  wavenumber  selection  problem  for  a  two- 
dimensional  geometry  (rolls),  and  in  the  case  where  the  Prandtl  number  can  be  considered 
infinite  (the  possibility  of  boundary-layer  instabilities  [20],  although  not  observed  in  our 
analysis,  is  also  briefly  discussed).  This  is  done  in  section  III. 3.  Buoyancy  effects  will  be 
neglected,  in  order  to  concentrate  on  the  effects  of  thermocapillarity.  Calculations  are 
achieved  for  a  semi-infinite  system  [14]  (i.e.  ignoring  the  presence  of  a  rigid  lower 
boundary,  and  thus  focusing  on  interfacial  short-wave  effects).  This  simplification  allows  us 
to  obtain  approximate  analytical  results  (section  III. 4)  about  the  convective  heat  transfer  far 
from  the  threshold,  and  about  other  relevant  quantities  such  as  interfacial  velocities  and 
surface  temperature  variations.  We  end  this  section  by  extending  some  of  these  results  to 
three-dimensional  disturbances,  and  reconsidering  the  problem  of  the  competition  between 
rolls  and  hexagonal  convective  structures. 


III.2.  Problem  formulation  -  weakly  nonlinear  results 

We  consider  a  semi-infinite  viscous  Boussinesquian  incompressible  fluid  in  contact  with  an 
inert  gas  phase.  The  interface  is  located  at  the  z=0  coordinate  plane  of  a  cartesian  reference 

frame  with  unit  vectors  1,.  ii=x,y,z),  and  is  assumed  undeformable  (this  will  allow 
obtention  of  analytical  results,  and  is  justified  since  interfacial  deformation  is  known  to 
primarily  affect  long-wave  modes  [21,22]).  The  fluid  is  located  in  the  domain  z<0,  and  a 
constant  heat  flux  is  injected  in  the  system  (a  constant  temperature  gradient  -13  is  maintained 
at  z^-oo).  All  equations  and  boundary  conditions  are  scaled  by  d  (an  arbitrary  length)  for 
length,  d'/K  for  time,  ^d  for  temperature  and  hkM^  for  pressure.  The  Marangoni  number 
Ma  =  -  Oj.fi d^/ IX  K  is  defined  with  respect  to  the  length  d,  instead  of  the  fluid  thickness  h 

(h/d-^oo).  Let  y  =  y^  +  W  be  the  fluid  velocity  (y^  is  the  horizontal  velocity),  T  the 
temperature  and  p  the  pressure  perturbations  with  respect  to  the  purely  conductive  (zero 
velocity)  reference  solution.  A  solution  vector  U  will  then  be  defined  by 


U{r=xl^+yly,Z,t) 


which  is  assumed  to  belong  to  a  certain  set,  say  E,  of  sufficiently  derivable  functions 


r 

w 

p 

T 


(3.1) 
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satisfying  the  boundary  conditions  of  the  problem  :  these  are 

V^,W,DT,p^0  forz^-o=  (3-2) 

W  =  DT  +  BiT  =  0  forz=0 

where  D  is  the  dimensionless  z-derivative  and  Bi  is  the  free  surface  Biot  number B/  =  ad l\ 
{a  is  the  free  surface  heat  transfer  coefficient  and  X  the  thermal  conductivity  of  the  fluid). 


The  system  of  partial  differential  equations  for  the  solution  vector  U  can  be  written  under  the 
general  operational  form  (see  eqs  2. 7-2. 9) 


X{U)  =  MaM(U)  +  Q(^)+NiU,U) 

ot 


(3.4) 


where  the  linear  part  ^(U)  is  given  by 


i£(C/)  = 


AV^-V^p 
AW -Dp 

AT+W 


the  "evolution  part"  is  defined  as 


Q{U)  = 


roi 

0 

0 

T 

0 


and  the  "constraint  part"  by 


M{U)  = 


Finally,  the  bilinear  form  N  is  expressed  as 


N{U, ,  U^) 


0 

0 

0 

0 


0 

0 

0 


(3.5) 


(3.6) 


(3.7) 


(3.8) 
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In  the  above  relations,  =  1  —  +  is  the  horizontal  gradient,  V  =  V  +  1  Z)  is  the  total 

dx  ^dy  '  ^ 

gradient,  and  A  =  is  the  Laplacian  operator. 

Apart  from  the  fact  that  pressure  is  not  directly  eliminated  from  the  equations,  the  set  (4)  of 
equations  together  with  boundary  conditions  (2-3)  is  equivalent  to  the  problem  formulation 
of  Scanlon  and  Segel  [14],  Note  also  that  we  have  included  the  Marangoni  condition  as  the 
last  component  of  (4),  which  has  already  been  shown  to  simplify  the  process  of  deriving 
amplitude  equations  [14,23-26], 

i)  Derivation  of  amplitude  equations 

We  first  decompose  U  into  Fourier  modes 

U(  r,z,t)  =  I  Uj{z,t)  exp(/  k.r)dk 

so  that  horizontal  Fourier  components  Uj  all  belong  to  E  (i.e.  fulfil  boundary  conditions  (2- 
3))  and  satisfy 

=  Ma  MjiUj)  .  Gji^)  +  j dV 


(3.9) 


(3.10) 


which  is  obtained  by  projecting  (4)  on  exp(-/^.  r)  and  by  replacing  by  ik  in  linear 

operators  (this  is  the  meaning  of  the  index  k).  The  bilinear  form  N  is  defined  in  a  similar 
way.  Each  Fourier  mode  is  further  decomposed  as 


U^zj)  =  A^t)  UKz)  +  U^{z,t) 


(3.11) 


where  Uj{z)  is  an  eigenvector  with  eigenvalue  o^ik  =  |  ^| )  of  the  linear  spectral  problem 

.  a,ej(U{)  (3.12) 

The  resolution  of  (12),  detailed  in  appendix  2,  shows  that  for  any  Ma>0,  0<k<Ma/2  Bi, 
an  isolated  eigenvalue  exists  (and  is  such  that  ai^+lc^>0).  This  eigenvalue  is  the  growth 

rate  of  the  corresponding  eigenmode  Uj{z),  appearing  in  eq.  (11).  For  every  value  of  Ma 
and  k,  there  also  exists  a  continuum  of  solutions  of  (12)  that  are  bounded  for  z-^-c»  (and 
which  correspond  to  eigenvalues  a<  =-k^.  This  infinite  set  of  solutions  could  eventually  be 

used  to  develop  the  remainder  term  Uf(z,t)  (the  superscript  D  stands  for  "damped"),  but  it 

turns  out  to  be  simpler  to  compute  Uf  directly,  by  a  method  explained  in  appendix  1.  As 
exchange  of  stability  holds  in  our  problem  [27] ,  eigenvalues  o)^  are  real  and  satisfy 


,,  2,  .2  o-^+2k^a-^+2kBia-\.i 

Ma  =  —  (  -  2ka  ^  + - )  '■ 

k 


Bi  +  \[a~¥k^ 


(3.13) 


The  neutral  stability  condition  is  found  by  the  limit  of  equation  (13)  for  a-^0  : 

Ma,  =  8k(k+Bi)  (3.14) 

which  is  the  asymptotic  form  (^->oo)  of  the  neutral  stability  condition  of  Pearson  [13],  as 
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seen  on  fig.  1.  A  relation  between  the  maximal  eigenvalue  and  Ma  may  be  found 

by  differentiating  (13)  at  constant  Ma,  setting  da/dk=0.  This  gives,  for  Bi-0  for  example: 

‘^max  =  =  ocMu  wUH  «  =  0.086  (Bi  =  0)  (3-15) 


Although  analytical  results  can  also  be  obtained  for  Bi^^O,  they  are  not  reproduced  here  for 
conciseness.  As  remarked  by  Scanlon  and  Segel  [14],  the  minimum  (critical)  value  of  Ma  for 
which  instability  occurs  is  zero,  due  to  the  absence  of  stabilisation  by  a  rigid  lower  boundary 
of  the  modes  with  increasingly  large  wavelength.  However,  due  to  their  large  inertia,  the 
growth  rate  of  these  modes  is  vanishingly  small  for  all  Marangoni  numbers.  This  is  depicted 
in  fig. 2  :  it  is  seen  that  modes  with  wavenumbers  between  0  and  k*=(Ma/8y^^  (for  Bi=0) 
are  unstable,  so  that  their  amplitude  Ayt)  in  the  decomposition  (11)  should  grow 
exponentially  in  time,  as  long  as  nonlinear  effects  can  be  neglected.  In  fact,  it  is  shown  in 
appendix  1  that  the  amplitudes  A^t)  obey  evolution  equations  of  the  form 


dAj 

dt 


These  equations  are  strictly  valid  near  threshold.  When  the  Marangoni  number  is  increased, 
higher  order  terms  should  be  included.  Equations  (16)  may  then  be  considered  as  resulting 
from  a  truncated  modified  Galerkin  scheme  [15,25].  Another  hypothesis  underlying  the 

derivation  of  (16)  is  that  the  dynamics  of  damped  modes  (i.e.  of  Uj)  is  determined  by  the 

evolution  of  the  "primary"  modes  Ayt)  U{  (this  amounts  to  neglect  time  derivatives  of 
damped  modes).  This  slaving  principle  [25,26,28],  strictly  valid  near  the  threshold,  is  here 
assumed  to  be  qualitatively  valid  in  the  strongly  non-linear  regime.  This  can  be  partly 
justified  by  the  fact  that  damped  modes  cannot  bifurcate  whatever  high  the  value 

of  Ma  (see  also  [33]). 


Despite  these  assumptions,  our  model  is  expected  to  reflect  physical  reality  even  far  from 
threshold,  provided  that  the  eigenmodes  Uj  are  used  for  the  Galerkin  basis,  rather  than  the 

neutral  stability  functions  U^.  In  order  to  illustrate  the  differences  between  these  different 
approaches,  we  now  turn  to  the  derivation  of  weakly  nonlinear  results  [14,  23-26],  for  which 
the  latter  option  is  sufficient. 


ii)  Weakly  nonlinear  results 


Making  use  of  Uyz,f)  =  Ayt)  Uj{z)  +  instead  of  eq.ll,  and  following  a  procedure 

similar  to  that  described  in  appendix  1 ,  we  are  left  with  amplitude  equations  identical  to 
eq.  (16),  although  with  different  coefficients.  It  is  obvious  that  coefficients  of  the  quadratic 
and  cubic  terms  do  not  depend  on  Ma.  It  can  also  be  shown  that  the  coefficient  of  the  linear 
term  is  the  first  term  of  a  Taylor  expansion  of  a^Ma)  around  Ma=Ma^,  i.e. 

0  _  kiMa-Ma,)  (3.17) 

~  A{2k+Bi) 

where  a  superscript  0  will  denote  a  value  of  a  coefficient  computed  by  using  neutral  stability 
functions. 
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Figure  3.2  :  The  growth  rate  a  as  a  function  of  the  wavenumber  k  for  different  Marangoni 
numbers  Ma,  and  for  Bi—0.  The  dotted  line  represents  the  locus  of  the  fastest  growing 
perturbations,  given  by  eq.  (15). 

Note  that  eq.  (17)  is  an  insufficient  approximation  for  low  Biot  numbers.  As  an  indication, 
the  mode  k=0  is  found  to  be  amplified  for  Bi-0  and  Ma>0  (which  differs  from  the  exact 
behaviour  of  the  growth  constant,  as  seen  in  fig. 2). 


From  equations  (9)  and  (11),  it  is  seen  that  a  roll  mode  with  wave  vector  k^  is  described  by 

A-j^  =  (k-  k^)  +a^(t)8(k+kQ),  where  8  is  the  Dirac-delta  function,  and  an  overbar 

denotes  the  complex  conjugate.  Substituting  into  (16)  leads  to  the  Ginzburg-Landau  equation 


^  +(2Z°„+z!„.)«, 


(3.18) 


where  Zi°  i  and  Z°,i,  stand  for  Z|^j^  and  Z_jj^  respectively. 


-  kg  k^ 


Defining  a  reduced  distance  to  the  threshold  by 

e  =  (Ma-MaAIMa,  (3.19) 

it  is  seen  that  at  e=0,  the  rest  state  aj=0  undergoes  a  pitchfork  bifurcation  to  the  steady 
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amplitude 


1 

I 


32{Bi  +  2>k^) 

{Bi  +k^)09Bi^+2Amk^  +353/to) 


obtained  after  evaluation  of  the  cubic  coefficients  (see  details  in  appendix  3). 


(3.20) 


Although  strictly  valid  near  the  threshold,  the  limitations  of  this  weakly  nonlinear  model  for 
large  e  are  well  known.  Consider  for  example  the  temperature  perturbation  averaged  in  the 
horizontal  direction  (i.e.  its  k=0  Fourier  component)  ; 


(r>  =  Tl,{z)  =  2\a,\^T^,{z)  =  6>(€) 


(3.21) 


where  7’o°(z)  is  the  only  non-zero  component  of  Uqx{z)  (appendix  3).  The  total  averaged 
temperature  profile  is  obtained  by  adding  the  reference  profile  -z  to  (21),  and  is  represented 
in  fig.  3,  for  several  values  of  e.  It  is  seen  that  in  a  region  of  depth  0(1 /k()  below  the 
interface,  the  temperature  profile  is  distorted  (and  somewhat  homogenised)  by  Marangoni 
convection.  It  is  also  seen  that  unrealistic  temperature  distributions  (strongly  negative  values 
of  the  mean  temperature,  leading  to  large  unrealistic  cold  spots  in  steady  regimes)  are 
obtained  for  e  superior  to  about  1.5.  Defining  A  as  the  bulk  temperature  decrease  with 
respect  to  its  value  in  the  conductive  rest  state,  we  may  compute  that 

.  32{Bi^3k,)(5Bi^lk,) 

|Mimroi(z)  - 


=  -2  I  a, 


k^  (39  BP + 24Sk^Bi  +  353  ko) 


(3.22) 


where  the  superscript  0  again  denotes  the  weakly  nonlinear  result.  The  result  (22)  of  course 
diverges  for  e-*oo. 


The  importance  of  obtaining  a  better  approximation  of  the  bulk  temperature  decrease  A  is 
justified  by  the  fact  that  it  can  be  considered  as  equivalent  to  the  classical  Nusselt  number 
Nu  (more  exactly  to  Nu-1,  which  is  also  quadratic  in  the  amplitudes).  Indeed,  for  systems 
in  which  the  temperature  difference  is  kept  constant  (such  as  Rayleigh-Benard  convection 
between  conducting  boundaries),  Nu  is  defined  as  the  dimensionless  ratio  of  the  total  to  the 
conductive  heat  flux,  and  therefore  is  a  measure  of  the  increase  of  the  heat  flux  due  to 
convection.  For  systems  where  the  heat  flux  is  kept  constant  (as  in  the  present  work),  the 
decrease  of  the  temperature  difference  between  bulk  and  interface  may  also  be  perceived  as 
an  increase  of  the  apparent  thermal  conductivity  of  the  system,  due  to  Marangoni-Benard 
convection.  In  the  next  sections,  it  is  shown  that  by  using  eigenvectors  (12)  instead  of  neutral 
stability  functions,  a  more  realistic  description  of  convective  fields  for  large  Ma  can  be 
obtained,  together  with  interesting  power  laws  for  the  variation  of  convective  quantities  in 
the  limit  Ma-^oo . 


III. 3.  Numerical  results  and  physical  interpretation 

In  this  section,  we  present  results  obtained  by  direct  numerical  integration  of  the 

set  (16),  for  a  two-dimensional  domain  of  lateral  length  L=2-K/ko  with  periodic  boundary 

conditions. 
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Figure  3.3  :  The  total  temperature  <  7Voj->  averaged  in  the  horizontal  plane  as  a  function 
of  the  vertical  coordinate  z,  as  computed  from  weakly  non-linear  results,  and  for  different 
values  of  e  =  (Ma-MaJ/Ma^.  The  Biot  number  is  Bi=0,  and  the  basic  wavenumber  is  ko=l. 
The  distortion  of  the  averaged  temperature  profile  in  the  convective  region  near  the  interface 
z=0  creates  a  decrease  Aq  of  the  bulk  temperature  with  respect  to  the  purely  conductive  value 
(dotted  line,  e=0).  Aq  is  here  defined  for  e= 0.5. 


The  amplitude  of  the  Fourier  modes  is  given  by 

=  Yl  aj,t)b{k-nk^) 

n=-N, 


(3.23) 


with  a,j(r)  =  and  N  sufficiently 

Substituting(23)  in  relation  (16)  leads  to 

dt 


large 


=  cr.„a„ 


.  a  a 

p,m  p  m-p 


to  ensure  numerical  convergence. 


^  ^  ^p,q,m^p^q^m-p-q 
P><I 


(3.24) 


where  a„,  stands  for  for  and  Z^  „,„,  for  Z^j^ 

a  function  of  the  Marangoni  number  Ma. 


which  are  calculated  as 


In  the  following,  we  will  take  advantage  of  the  fact  that  the  length  scale  d  of  the  problem  is 
still  arbitrary.  We  may  thus  choose  ko=l,  which  means  that  the  dimensional  length  of  the 
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periodic  box  is  lird.  From  eq.  (14),  the  critical  Marangoni  number  is  given  by 
Ma^=Mai^^j—8  (1+Bi). 

The  system  of  equations  (24)  has  been  integrated  for  a  wide  range  of  Marangoni  and  Biot 
numbers.  Despite  the  large  number  of  unstable  modes  in  some  cases  (increasing  with  Ma), 
and  the  presence  of  resonant  quadratic  terms  (which  are  generally  responsible  for  complicated 
phase  coupling  effects  [29,30]),  the  long-term  behaviour  appears  to  be  surprisingly  simple: 
a  steady  state  is  always  reached,  which  is  strongly  dominated  by  the  fundamental  mode  «=i. 
Since  the  number  of  modes  N  needed  to  ensure  convergence  is  increasing  with  Ma  (see 
fig. 2),  computer  resources  limited  our  investigations  to  Marangoni  numbers  of  about 
Ma=500  (for  Bi=0,  and  N=20). 


We  have  also  considered  a  simplified  version  of  the  system  (24),  which  allows  us  to  simulate 
the  evolution  of  a  larger  number  of  amplitudes.  This  model  is  obtained  by  setting  all  cubic 
coefficients  Z„„,„  with p7^m  equal  to  zero.  We  then  obtain 

aa..  ^ 


—!l  =  [a 
dt 


y  s  \a  \^]a  +  y  Z 

/  ^  m,q  I  ^  I  -*  ni  /  ^ 


a  a 

p,m  p  m-p 


(3.25) 


with  S,„  ^=-2  Z„,  ^  From  eq.  (25),  this  quantity  is  seen  to  represent  the  strength  with  which 
the  presence  of  the  mode  q  lowers  the  effective  growth  rate  of  the  mode  m.  The  physical 
mechanism  responsible  for  this  stabilising  effect  consists  in  the  distortion  of  the  mean 
temperature  profile  by  convection  (see  fig. 3),  which  lowers  the  destabilising  temperature 
gradient.  A  comparison  of  the  time  evolution  and  of  the  steady  state  values  predicted  by  the 
full  system  (24)  and  the  reduced  set  (25)  reveals  that  the  results  differ  only  slightly  (by  less 
than  10%  on  the  value  of  typical  convective  quantities  at  steady  state,  as  shown  in  fig.  6). 
In  view  of  this  rather  good  concordance,  the  mean  temperature  profile  distortion  by 
convection  may  be  considered  as  a  dominating  effect  in  the  nonlinear  competition  between 
unstable  modes.  Implications  for  wavelength  selection  between  fastest  growing  and  critical 
modes  are  discussed  later  on  in  this  section. 


Since  it  is  legitimate  to  admit  that  the  simplified  system  (25),  that  can  be  considered  as  a 
"mean-field"  [32]  version  of  the  problem,  is  useful  for  simulating  the  interactions  of  a  larger 
number  of  modes  (up  to  N=75),  larger  Marangoni  number  values  can  also  be  investigated. 
Again,  even  for  Marangoni  numbers  as  large  as  4000  (for  Bi—0,  i.e.  e  —  500),  the  long-term 
behaviour  is  not  modified,  independently  of  the  initial  conditions  (here  selected  as 
a  numerical  "white  noise",  i.e.  randomly  chosen  complex  amplitudes  of  magnitude  to 
7(7^)  :  the  final  state  is  still  steady  and  dominated  by  the  fundamental  mode. 

A  sequence  of  a  typical  transient  simulation  is  represented  on  fig. 4.  For  sufficiently  small 
initial  perturbations,  a  convective  structure  dominated  by  the  fastest  linearly  growing  mode 
(the  mode  closest  to  is  observed  after  a  relatively  short  time.  This  is  the  case  as  long 
as  non-linear  effects  can  be  neglected.  At  higher  time  intervals,  this  ^,,,^,-structure  is 
progressively  replaced  by  larger  and  larger  wavelength  structures,  via  a  complex  process  of 
coalescence  of  neighbouring  convective  cells.  This  evolution  finally  tends  to  the  steady  state 
with  2  convective  cells  (1  period)  occupying  the  entire  domain,  as  expected  (fundamental 
mode).  The  properties  of  this  steady  state  will  be  investigated  in  the  next  section. 
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Figure  3.4  :  The  evolution  of  the  streamfiinction  pattern  for  Ma=3000  (e=374),  Bi=0, 
N=75.  Solid  closed  curves  correspond  to  clockwise  motion.  The  initial  condition  was 
selected  as  a  random  noise  of  amplitude  The  streamfiinction  is  rescaled  at  each  snapshot 
(the  reduced  time  is  indicated).  The  fastest  growing  perturbation  is  dominating  for  times 
t<0.03.  A  continuous  growth  of  the  mean  wavelength  of  the  pattern  is  observed,  the  later 
stages  of  which  tend  to  a  steady  state  with  two  convective  cells  (1  period)  in  the  simulation 
domain,  after  a  time  of  order  unity  (in  units  where  L  is  the  horizontal  period). 
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It  is  interesting  to  compare  fig. 4  with  the  experimental  observations  made  by  Linde  [5],  in 
mass  transfer  systems.  Note  that  these  systems  actually  correspond  to  high  values  of  Bi,  since 
the  diffusion  coefficients  are  generally  much  larger  in  the  gas  than  in  the  liquid  phase. 
However,  our  simulations  were  not  found  strongly  dependent  on  Bi  for  large  Ma.  Linde 
interpreted  the  observed  growth  of  the  mean  wavelength  of  the  convective  pattern  as  an  effect 
due  to  the  non-stationary  mass  transfer  occurring  in  his  experiments.  Indeed,  after  that  the 
experiment  is  started  by  putting  in  contact  a  gas  phase  containing  a  surface  active  solute  with 
the  liquid  phase,  the  diffusion  of  this  solute  through  the  interface  creates  a  growing  diffusive 
boundary  layer,  which  induces  convective  motions  in  the  liquid,  with  a  wavelength  that  scales 
with  the  thiekness  of  this  boundary  layer. 

Since  a  natural  length  scale  such  as  this  boundary  layer  thickness  is  absent  in  our 
formulation,  the  wavelength  selection  observed  in  fig. 4  has  to  be  intrinsically  related  to  the 
non-linear  mechanism  of  heat  (or  mass)  convective  transport.  This  effect  was  indeed  shown 
(see  fig. 3)  to  create  an  homogenisation  of  the  temperature  (or  concentration)  profile  in  a 
convective  region  located  below  the  interface.  This  is  also  apparent  in  fig. 5,  which  represents 
the  temperature  profile  averaged  along  the  horizontal  direction  corresponding  to  the  evolution 
depicted  in  fig. 4.  It  is  seen  that  the  temperature  uniformisation  due  to  convection  is  more 
important  at  large  times,  when  the  penetration  depth  is  large.  The  growth  mechanism  can  be 
explained  by  the  following  considerations  :  suppose  that  at  one  particular  instant,  the 
convective  structure  has  a  given  mean  wavelength  X.  Since  the  convective  cells  have  to 
preserve  a  certain  height/width  ratio,  temperature  is  practically  homogenised  in  a  region  of 
depth  X  below  the  interface.  Modes  with  wavelengths  smaller  than  X  may  be  considered  as 
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Figure  3.5  :  The  evolution  of  the  total  temperature  <  Tjot>  averaged  in  the  horizontal  plane 
as  a  function  of  the  vertical  coordinate  z-  The  Biot  number  is  Bi=0.  Several  times  are 
considered,  which  correspond  to  the  simulation  depicted  in  fig. 4.  An  homogenisation  of  the 
temperature  in  a  domain  whose  depth  is  growing  with  time  occurs,  which  results  in  a  growth 
of  the  bulk  temperature  decrease  A  with  time  (see  insert). 


stable,  since  they  see  a  nearly  isothermal  environment.  At  the  contrary,  modes  with 
wavelengths  larger  than  X  can  penetrate  deeply  enough  into  the  bulk  of  the  liquid  and  bring 
hot  fluid  from  the  still  conductive  zone,  to  the  interface.  The  effective  growth  rate  of  these 
modes  remains  nearly  unchanged  by  the  convective  structure,  so  that  these  modes  continue 
to  grow  (but  slower  and  slower  due  to  their  growing  inertia),  and  tend  to  replace  smaller 
wavelength  structures.  The  pattern  wavelength  X  may  thus  be  expected  to  grow  indefinitely, 
at  least  in  an  infinite  system.  In  real  experiments,  the  final  wavelength  will  probably  be 
determined  by  the  actual  depth  of  the  experimental  container  (thus  near  the  critical 
wavelength),  indicating  that  an  intrinsic  wavelength  is  likely  to  be  inexistent  for  the  pure 
Marangoni-Benard  problem.  Note  that  on  the  point  of  view  of  wavelength  selection,  the 
evolution  described  above  presents  some  similarities  with  coarsening  processes  observed 
during  spinodal  decomposition  phenomena  in  binary  mixtures  [36]. 

A  last  remark  about  fig. 5  concerns  the  temperature  profile  near  the  interface.  Since  the 
vertical  velocity  is  vanishing  at  the  interface,  some  kind  of  thermal  boundary  layer  is  created 
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there,  in  which  the  temperature  gradient  quickly  recovers  its  bulk  value.  In  Rayleigh-Benard 
convection,  boundary  layer  effects  are  known  to  play  a  decisive  role  in  the  mechanisms  of 
transition  to  turbulence  (especially  for  high  Prandtl  number  fluids  [20]).  However,  despite 
the  very  high  values  of  the  Marangoni  number,  boundary  layer  instabilities  have  not  been 
observed  in  our  simulations,  probably  due  to  the  different  nature  of  these  boundary  layers 
(in  particular  the  absence  of  no-slip  condition  for  Marangoni-Benard  problems).  Furthermore, 
it  cannot  be  rejected  that  this  kind  of  phenomena  could  appear  for  larger  driving  forces  than 
those  investigated  in  this  work  (up  to  Ma=4000).  Finally,  let  us  mention  that  a  direct 
comparison  of  the  results  obtained  from  the  present  model  (amplitude  equations  limited  to 
third  order)  with  a  finite  difference  resolution  of  the  governing  equations  is  in  progress,  and 
will  be  reported  elsewhere.  Preliminary  results  exhibit  a  satisfactory  agreement  conce^ng 
the  qualitative  evolution  of  the  system  (i.e.  the  growth  of  the  mean  wavelength  up  to  the  final 
steady  state  with  the  largest  wavelength).  This  confirms  that  the  most  important  ingredient 
responsible  for  this  process  is  indeed  the  mean  temperature  profile  distortion  by  convection. 
This  in  turn  indicates  that  mean-field  approximations  [32],  neglecting  all  nonlinear  effects 
except  the  change  of  the  mean  temperature  profile  owing  to  the  convective  heat  transport, 
can  lead  to  satisfactory  approximations  of  highly  supercritical  behaviours. 

III.4.  Analysis  of  steady  states 

i)  Bifurcation  of  rolls 


Since  the  steady  state  reached  by  both  full  (24)  and  reduced  (25)  models  is  strongly 
dominated  by  the  fundamental  mode  n=l,  we  seek  for  an  approximate  solution  by  setting 
to  zero  all  harmonics  with  1.  The  set  (25)  then  reduces  to  the  single  equation 


da^ 

It 


[(7-5i,  |aJ2]ai 


(3.26) 


describing  a  pitchfork  bifurcation  similar  to  eq.  (18)  but  where  the  coefficients  are  now 

computed  from  the  eigenfunctions  17/,  and  thus  depend  on  Ma.  After  computation  of  these 
coefficients,  the  steady  convective  solution  of  (26)  is  finally  found  as 

|^j2  ^  ^  _ (Ma-Ma^)<7^(3+v/l+ff  f _  ^7) 

(1  +5/)^(512(A/<3+o^)  +  A/(3(3cr-8)(3+y^l  +0 )^) 


where  the  growth  rate  a  is  solution  of  the  dispersion  relation  (13),  and  thus  also  depends  on 
Ma.  According  to  eq.  (26),  the  solution  (27)  is  stable  provided  a>0,  which  is  equivalent  to 
Ma>Mac=8(l+Bi)  (we  have  s&t  ko=l). 


The  decrease  A  of  the  bulk  temperature  due  to  Marangoni  convection,  as  represented  in  figs 
3  and  5,  is  expressed  by 

8(Ma-Ma^)(3+VW 8(Ma+o^)  ^ 

^  ^  _ Majl+sJUp  f  (3.28) 

5ll{Ma+o^)  +  Mfl(3(T-8)(3+\/l+a )^ 

This  expression  is  represented  in  fig. 6,  together  with  results  obtained  from  the  integration 
of  the  full  system  (24)  and  of  the  reduced  system  (25).  Another  result  found  in  fig. 6  is  the 
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expression  (22),  which  reduces  to  A°=672e/353  for  Bi=0  and  ko=l.  Clearly,  this  result  is 
only  valid  near  the  threshold.  At  the  contrary,  the  expression  (28)  for  A  leads  to  the  result 
A=14/5  e  near  the  threshold  (which  is  overestimated,  due  to  the  negligence  of  the  stabilising 
coefficient  see  eq.l8).  Nevertheless,  eq.  (28)  appears  to  be  a  better  approximation  of 
the  bulk  temperature  decrease  for  large  Marangoni  numbers  (because  Z_jjj  becomes  negligible 
compared  to  The  corresponding  mean  temperature  profile  can  also  be  shown  to  be  more 
realistic,  since  it  does  not  exhibit  cold  spots  as  those  appearing  in  fig. 3,  but  is  rather  close 
to  fig. 5.  The  behaviour  at  an  infinite  Marangoni  number  is  given  at  the  end  of  this  section. 


Figure  3.6  :  The  bulk  temperature  decrease  A  as  a  function  of  the  distance  to  the  threshold 
e  =  (Ma-MaJ/Ma^  for  Bi=0.  Thick  full  curve  :  results  of  the  numerical  integration  of  the  full 
system  (25).  Thick  dotted  curve  :  numerical  integration  of  the  "mean-field"  system  (24).  Thin 
full  line  :  the  weakly  nonlinear  result  A^=672e/353.  Thin  dotted  curve  :  the  analytical  result 
for  A  given  as  eq.  28  of  the  text.  The  insert  represents  a  zoom  of  a  region  near  the  origin. 

ii)  Competition  between  hexagons  and  rolls 

In  view  of  the  good  agreement  between  the  analytical  result  (28)  and  the  results  of  the 
numerical  integration  of  (24),  we  shall  reexamine  the  problem  of  the  competition  between 
three  sets  of  rolls  forming  angles  of  60°  with  each  other. 

We  thus  consider 

A-j;  =  a^(f)8(k-k^)  +  a^{f)b{k-k^  +  a^(t)8{k-k^+ 
a^{t)b{k+k^  +  a^(t)b(k+k^)  +  a^(t)bik+k^) 


(3.29) 


where  the  ordering  of  unit  vectors  ,  |/|  =  1,2,3  is  defined  by  fig.7. 


Figure  3.7  :  Definition  of  the  basic  wavevectors  for  the  study  of  the  competition  between 
rolls  and  hexagons  ( |  ^,|  =1). 


From  eq. 


where 


(16),  the  corresponding  amplitude  equations  are 


dt 

=  oa^ 

+ 

6(32 

a,  - 

0!,  |a,|2 

+  Q!2( 

l«2 

1^-! 

i«3r)' 

da^ 

dt 

8a^ 

^3  - 

a, \a2V 

+ 

|a, 

l«3l^)' 

^2 

(3.30) 

da^ 

dt 

+ 

6  a, 

«2  - 

+  a2( 

|2+ 1 

6 

=  2Z,2 

(3.31) 

a, 

=  “  (2^1,1, 1  + 

^-l.U 

) 

(3.32) 

a. 

2 

-2(Z 

1,1,1  ^2,1 

■2,1,1 

) 

(3.33) 

and  where  Z^  ^  ,.  stands  for  Zjjj  (symmetry  considerations  have  been  used  to  minimize  the 

number  of  coefficients  to  be  calculated).  The  discussion  of  the  gradient  system  (30)  is  well- 
known  [23-26,37]  :  writing  a„=r„  exp[i<pj,  we  obtain  the  equation 

d<p/dt  =  -8  sin^j  (r^r^  +rlrl  +r1rl)lr^r^r^  (with  0)  for  <p=(p2-<Pr^3,  showing  that  <p=0 
and  <p  =  'K  are  the  only  possible  stationary  values  of  cp. 
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Figure  3.8  :  The  coefficients  b/(l  +Bi)  (upper  graph),  a/(i  (lower  graph,  solid  curves) 
and  a2/(l  -VBif  (lower  graph,  dashed  curves)  as  a  function  of  the  distance  to  the  threshold 
e  =  (Ma-MaJBi))/MaJBi),  for  various  Biot  numbers  Bi  (indicated  on  each  curve). 
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Then,  it  is  found  that  qualitatively  different  fixed  points  of  eqs  30  are  (when  o:=o:;+2q;2>  0): 

-  rest  solution  ;  rj  =  r2  =  =  0  (3.34) 

-  roll  solutions  :  Tj  =  =  0.  =  {oloLiy^  (3.35) 

-  up-hexagons  solutions  ;  (p  =  0  ,  =  r2  =  •\-4aaf'^)/2oi.  (3.36) 

-  down-hexagons  solutions  :  <p  =  ir  ,  =  r2  =  =  (-h  +  (b^ +4aay'^)/2a  (3.37) 

The  analytical  form  of  the  coefficients  b,  a j,  and  a2  (depending  onMa)  is  not  written  down 
for  conciseness.  Rather,  fig.  8  presents  their  variation  with  the  distance  to  the  threshold  e  for 
various  Biot  numbers.  Bifurcation  diagrams  are  represented  on  figs  9  and  10. 


Figure  3.9  :  The  bulk  temperature  decrease  A  as  a  function  of  the  distance  to  the  threshold 
e  for  Bi=0.  R  :  Rolls,  U-H  :  Up-hexagons,  D-H  :  Down-hexagons.  Solid  curves  represent 
stable  states,  while  dotted  curves  represent  unstable  states.  The  thin  solid  line  represents  the 
analytical  result  given  by  eq.  (28)  of  the  text  for  the  bulk  temperature  decrease  of  rolls. 


Figure  9  represents  the  bulk  temperature  decrease  A  as  a  function  of  the  distance  to  the 
threshold  e,  for  solutions  (35)  to  (37).  As  expected,  up-hexagons  (upflows  at  the  centre  of 
the  hexagons)  are  the  only  stable  solutions  just  above  the  threshold,  and  rolls  become  stable 
only  at  large  €.  Down-hexagons  (downflow  at  the  centre)  are  always  unstable,  because  b>0 
(which  is  different  from  the  case  of  very  low-Prandtl  number  fluids  [26]).  Note  that,  although 
not  apparent  on  the  figure,  the  first  bifurcation  to  up-hexagons  at  e=0  is  slightly  hysteretic: 


47 


the  depth  of  this  subcritical  region  is  3.37o  (for  Bi=0)  in  our  model,  slightly  larger  than  the 
2.3%  value  of  Scanlon  and  Segel  [14].  This  is  due  to  the  fact  that  we  have  neglected  the 
stabilising  action  of  "secondary"  modes,  i.e.  those  generated  by  quadratic  interaction  of  the 
"primary"  modes.  This  is  done  since  it  is  natural  to  expect  that  these  modes  become 
unimportant  for  large  Ma,  as  observed  in  the  case  of  two-dimensional  simulations, 
characterised  by  a  strong  domination  of  the  fundamental  mode  (see  section  III. 3).  When 
the  amplitudes  of  these  harmonics  are  included  in  eq.  (29),  and  finally  eliminated  using 
adiabatic  slaving  [25-26,  28],  the  2.3%>  value  of  the  hysteresis  is  recovered.  It  is  also 
possible  to  recover  the  result  0.56%  of  Bragard  and  Lebon  [24],  in  the  case  of  a  layer  of 
finite  depth  (the  calculation  of  coefficients  is  then  fully  numerical).  It  is  also  apparent  that 
the  depth  of  the  subcritical  region  is  increasing  with  the  Biot  number. 

An  interesting  result  of  the  present  analysis  is  that  at  large  Marangoni  numbers,  the  stability 
properties  are  not  qualitatively  modified  with  respect  to  the  results  of  Scanlon  and  Segel  [14]: 
these  authors  predict  that  rolls  should  become  stable  above  a  value  ej=64  (our  value  is 
ej=8.6)  of  the  constraint,  while  up-hexagons  should  become  unstable  above  €2=196  (our 
value  is  €2=37).  A  bistability  region  (leading  to  hysteresis  effects  between  rolls  and 
hexagons)  thus  exists  between  and  €2.  This  qualitative  concordance  reinforces  the  idea  that 
this  hysteresis  region  could  be  a  physical  reality,  although  the  domain  of  validity  of  the 
amplitude  equations  is  not  guaranteed  for  such  large  values  of  e.  Lastly,  we  have  represented 
the  maximal  surface  velocity  for  the  bifurcating  solutions  (35)  to  (37)  in  fig.  10,  showing  that 
this  quantity  is  not  strongly  dependent  on  the  particular  planform  selected. 
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Figure  3.10  :  The  maximal  surface  velocity  as  a  function  of  the  distance  to  the  threshold 
e  for  Bi=0.  R  :  Rolls,  U-H  :  Up-hexagons,  D-H  :  Down-hexagons.  Solid  curves  represent 
stable  states,  while  dotted  curves  correspond  to  unstable  states.  The  thin  solid  line  represents 
the  analytical  result  for  the  maximal  surface  velocity  of  rolls. 


iii)  Asymptotic  behaviours  for  Ma-»oo 

It  results  from  the  examination  of  the  previous  figures  that  the  asymptotic  behaviour  of  the 
relevant  convective  quantities  for  Ma-^oo  obey  to  different  power  laws  than  those  generally 
derived  near  the  threshold  by  using  classical  perturbation  methods.  In  particular,  a  saturation 
is  observed  for  the  bulk  temperature  decrease  A  (which,  as  mentioned  earlier,  can  be 
considered  as  equivalent  to  the  Nusselt  number).  The  purpose  of  the  following  calculations 
is  to  derive  asymptotic  results  for  A,  for  the  maximal  surface  velocity  V^,  as  well  as  for  the 
surface  temperature  deviation  which  we  define  as  the  difference  between  maxima  and 
minima  of  temperamre  on  the  free  surface. 

From  the  dispersion  relation  (13),  it  is  straightforward  that 


kMa 


for  Ma  -»•  CX5 


(3.38) 
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(3.39) 


The  asymptotic  value  of  the  amplitude  of  rolls  is  derived  from  eq.  (27)  : 

1 

^ _ 

111  11 

2^  3^  (1+5/) 

and  it  is  calculated  from  eq.  (28)  that  the  saturation  value  of  the  bulk  temperature  decrease 
for  a  roll  structure  is 

^Rous  ^  8  ^  2.66  (3.40) 

3 

independently  of  the  value  of  the  Biot  number.  This  result  is  confirmed  by  all  the  curves  of 
fig. 6. 

It  is  readily  computed  that  the  maximal  surface  velocity  is 

=  3.61  Ma^  (3.41) 

2^3" 

while  the  amplitude  of  the  surface  temperature  variations  is  given  by 

_2 

->  =  14.66M«'^  (3.42) 

2^3^ 

The  corresponding  expressions  for  hexagons  (although  unstable  with  respect  to  roll 
dismrbances  for  Ma-^oo)  can  also  be  derived  analytically.  We  obtain 

(3.43) 

9 

3(1+5/)^  (3.44) 

CY  -  rv,  ^  (104733031  -60445052^/1)  ^  11.56(1+5/)^  (3-45) 

243(7+4v^) 

It  follows  that 

^Hexagons  -»  2.08  (3.46) 

for  the  bulk  temperature  decrease  of  hexagons  (both  up  and  down  hexagons  lead  to  the  same 
value  for  Ma-^oo,  as  seen  from  eqs  (36)  and  (37),  and  in  fig. 9).  This  value  is  inferior  to  the 
value  8/3  of  rolls.  However,  this  should  not  be  taken  as  a  rigorous  justification  for  the 
instability  of  hexagons,  since  it  well  known  that  the  principle  of  maximisation  of  the 
convective  heat  transport,  originally  proposed  by  Malkus  [31],  does  not  lead  to  stability 
predictions  that  are  generally  valid  [32] . 
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Finally,  the  maximal  surface  velocity  tends  to 

y^asor^  -  3.29 

and  the  amplimde  of  the  surface  temperature  deviations  is 

^  16,81 


(3.47) 

(3.48) 


To  check  the  assumptions  used  in  our  analysis,  it  should  be  interesting  to  compare  the  above 
results  with  a  full  numerical  integration  of  the  problem.  Preliminary  finite  differences 
simulations  indicate  a  slow  growth  of  the  bulk  temperature  decrease,  as  well  as  an  increase 
of  the  velocities  (although  with  an  exponent  larger  than  1/3)  coexistent  with  a  slow  decrease 
of  the  surface  temperature  variations.  This  allows  us  to  place  some  confidence  in  our 
analysis.  In  particular,  the  exponents  1/3  and  -2/3  may  be  considered  as  first  approximations, 
that  could  be  refined  by  deriving  higher  order  contributions  to  the  amplitude  equations.  It 
appears  also  that  the  present  asymptotic  analysis  leads  to  results  that  are  difficult  to  check 
by  finite  differences  simulations.  Indeed,  since  convergence  of  results  may  only  be  expected 
for  very  high  values  of  the  Marangoni  number,  numerical  difficulties  are  encountered,  mainly 
due  to  the  presence  of  very  steep  surface  temperature  gradients  at  the  cold  points,  where  the 
fluid  moves  downwards. 


We  conclude  this  section  by  remarking  about  an  important  mathematical  aspect  of  the  pure 
Marangoni-Benard  instability.  A  particular  feature  of  this  problem  is  that  the  neutral  stability 
condition  provides  Ma  as  a  single-valued  function  of  the  wavenumber  k.  This  means  that 
above  the  corresponding  critical  value,  one  and  only  one  eigenmode  is  unstable,  whatever 
large  is  the  value  of  Ma.  This  has  to  be  contrasted  [33]  with  Rayleigh-Benard  instabilities, 
where  n  eigenmodes  are  linearly  unstable  above  the  value  Ra„-^+n^-i{^)^/k^  of  the  Rayleigh 
number  (case  of  pure  Rayleigh-Benard  convection  between  stress-free  boundaries).  Clearly, 
the  above  analysis  would  require  non-trivial  modifications  to  account  for  interactions  between 
these  unstable  vertical  modes.  Physically,  this  difference  between  Marangoni-Benard  and 
Rayleigh-Benard  problems  is  certainly  related  to  the  different  natures  of  the  surface  and  the 
bulk  forces.  This  could  also  explain  why  neither  boundary  layer  instabilities,  nor  subsequent 
transitions  to  turbulence  have  been  observed  in  our  model,  in  the  range  of  Marangoni 
numbers  investigated. 
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III.5.  Conclusions  of  section  III 


Pure  thermocapillary  instability  in  layers  with  non-deformable  interface  and  infinite  Prandtl 
number  has  been  studied.  The  model  is  based  on  the  assumptions  that  the  dynamics  is 
determined  by  the  interactions  of  the  unstable  eigenmodes  of  the  linear  stability  problem,  and 
that  the  evolution  equations  describing  their  interactions  can  be  limited  to  third  order  in  their 
amplitudes.  Although  strictly  valid  near  the  threshold,  this  model  has  to  be  considered  as  an 
approximation  (similar  to  a  truncated  modified  Galerkin  method)  far  from  it.  Attention  has 
been  restricted  to  short-wave  effects,  so  that  the  layer  has  been  assumed  infinitely  deep. 

The  first  part  of  the  analysis  has  focused  on  the  two-dimensional  wavelength  selection 
problem  at  moderately  large  Marangoni  numbers,  for  which  the  wavelength  of  the  fastest 
growing  disturbance  is  much  smaller  than  the  critical  wavelength.  The  transient  numerical 
integration  of  third  order  amplitude  equations  has  shown  that,  although  a  structure  dominated 
by  the  fastest  growing  disturbance  can  appear  at  a  given  instant,  it  is  progressively  replaced 
by  larger  wavelength  structures.  A  steady  state  is  always  reached,  whose  wavelength  is  equal 
to  the  size  of  the  (periodic)  simulation  domain.  This  is  confirmed  by  a  finite  difference 
integration  of  the  problem.  Although  the  fastest  growing  wavelength  is  a  finite  quantity,  it 
appears  that  the  Marangoni-Benard  instability  cannot  induce  stable  localised  structures  with 
intrinsic  wavelengths  independent  of  the  experiment  dimensions  (at  least  up  to  Ma=4000). 
The  physical  mechanism  responsible  for  this  convective  cells  growth  process  is  shown  to  be 
related  to  the  nonlinear  convective  heat  transport,  which  creates  a  distortion  of  the 
temperature  field  in  the  region  located  near  the  free  surface.  The  distortion  of  the 
horizontally  averaged  temperature  profile  (the  cause  of  the  instability)  has  a  stabilising  effect 
on  short  wavelength  modes,  but  leaves  the  growth  rate  of  large  wavelength  structures 
relatively  unchanged.  This  produces  a  growth  of  the  pattern  wavelength,  which  presents  some 
resemblance  with  experimental  results  of  Linde  [5],  although  we  have  not  studied  this  point 
in  details,  since  transient  evolution  of  the  diffusive  (boundary  layer)  profile  is  not  taken  into 
account  in  our  analysis. 

Properties  of  the  steady  states  observed  in  the  convective  system  have  been  investigated  for 
both  rolls  and  hexagonal  structures.  Contrary  to  existing  weakly  nonlinear  theories,  our 
method  (differing  by  the  use  of  eigenfunctions,  rather  than  neutral  stability  functions)  appears 
to  lead  to  physically  realistic  results,  at  least  qualitative,  for  very  large  Marangoni  numbers. 
This  is  conjectured  from  examination  of  the  behaviour  of  some  relevant  convective  quantities. 
The  decrease  of  the  bulk  temperature  due  to  Marangoni  convection  has  been  found  to  present 
a  saturation  when  the  Marangoni  number  is  increased  (while  the  weakly  nonlinear  result  is 
a  linear  growth).  The  velocities  are  found  to  grow  as  as  Ma-^oo,  while  the  surface 
temperature  variations  decrease  as  However,  due  to  the  assumptions  underlying  our 

model,  these  behaviours  should  be  considered  as  first  approximations,  that  could  be  refined 
by  including  higher  order  interactions  in  the  amplitude  equations.  The  analysis  of  the 
competition  between  rolls  and  hexagons  confirms  earlier  results  [14,  23-26],  quantitatively 
(near  threshold)  and  qualitatively  (far  from  threshold).  At  small  Marangoni  numbers,  up- 
hexagons  are  the  only  stable  solutions  (the  first  bifurcation  is  slightly  hysteretic),  while  rolls 
are  the  only  stable  solutions  at  very  large  Marangoni  numbers.  Down-hexagons  are  always 
unstable  (Pr-^oo).  The  transition  between  up-hexagons  and  rolls  is  found  to  be  hysteretic,  as 
in  other  problems  (f.ex.  [34]).  However,  this  transition  has  been  found  to  occur  for  lower 
Marangoni  numbers  than  in  [14],  and  with  a  smaller  hysteresis  loop. 
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Appendix  1  :  Derivation  of  amplitude  equations 


Using  eq.  (11),  we  may  write  eq.  (10)  as  a  differential  problem  for 

dAr. 


'  -4?  (f'f .  Kt)  *  .?  ( .  Vi )) 


(A.  I) 


where  the  time-derivative  of  Uf  has  been  cancelled,  as  a  result  of  our  assumption  to  neglect 
the  own  dynamics  of  the  damped  modes  (slaving  principle).  Note  that  because  of  eq.  (11), 

the  boundary  conditions  can  only  be  satisfied  if  Uf  belongs  to  E.  A  very  rough  (and 
insufficient)  model  could  be  obtained  at  this  stage  by  projection  of  eq.  (A.l)  on  some 

functions  (generally  the  adjoint  functions  [35]),  and  assuming  Uf  =  0.  We  would  then 
obtain  an  equation  of  the  form  (16),  but  without  stabilising  cubic  terms.  Rather,  we  will  try 

to  solve  (A.l)  for  Uj  .  This  can  be  done  only  if  (A.l)  is  compatible,  which  is  not  the  case 
at  Ma—Ma^.  (the  kernel  of  the  left  hand  side  operator  is  then  not  empty),  except  if  the  second 
member  is  orthogonal  to  the  solution  of  the  adjoint  neutral  stability  problem  (Fredholm’s 
condition).  This  leads  to  the  amplitude  equations 


dAj 

dt 


I  ^k'  k^k'^k-k' 


I 


(A.2) 


where  eq.  (12)  has  been  used,  and  it  has  been  anticipated  that  the  velocity  (and  pressure) 
components  of  Uj  are  zero.  This  will  become  apparent  later,  and  is  a  consequence  of  the 

linearity  of  the  equations  of  motion.  In  (A.2),  (Vj  ,  denotes  the  projection  on  the  adjoint 

neutral  stability  solution  (derived  in  appendix  2),  =  (v^  is  the  normalisation 

factor,  and  the  quadratic  coefficients  are  given  by 

{vi  .'^TJ-viUi.Ulr)) 


^TI  - 


(A.3) 


Now,  by  inserting  the  projected  part  (A.2)  into  the  complete  equation  (A.l),  we  obtain 


k,k~k' 


(A.4) 


where  the  "Non-Resonant"  part  of  a  term  X  is  defined  by  X'^  =  X  - 

such  that  =  0  for  every  X.  This  ensures  that  (A.4)  is  compatible  for  every  Ma. 
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Now,  the  form  of  equation  (A.4)  suggests  an  iterative  series  solution,  starting  with' 

(A.5) 

in  which  U~,  is  obtained  from 


(A.6) 

Substituting  this  result  in  the  r-h-s  of  (A.4)  leads  to  higher  order  corrections  (cubic  terms  for 
Uj  ,  generating  quartic  and  higher  order  terms  in  (A. 2),  which  are  not  considered  in  our 

model).  We  limit  the  calculation  of  Uj  to  (A.5)  and  (A.6),  such  that  the  compatibility 
equations  (A. 2)  reduce  to  the  amplitude  equations  (16),  with  quadratic  coefficients  given  by 
(A. 3)  and  cubic  coefficients  given  by 

,  (Vr'.Al'.wfC'f.C'w.I-))  (A.7) 

~  ~ - - - 


Appendix  2  :  Solution  of  the  linear  problem 

Starting  from  eq.  (12),  the  neutral  stability  problem  {a^=0)  can  be  written  as 


SM)  =  X-M)  = 


{D^-k^)Wl-Dpl 

DWt+ik^j 


DV^j  +  ikMajl 


(A.8) 


The  solution  of  this  problem  which  belongs  to  E  (i.e.  which  satisfies  boundary  conditions  (2) 
and  (3))  reads 


Ir-] 

^  rk 

-Aik{k+Bi){\+kz) 

< 

= 

-Ak'^{k+Bi)z 

pI 

-Sk^k+Bi) 

/y^0 

1  -  (k+Bi)z  +  k(k+Bi)z^ 

where  k  =  |  )t|  ,  and  the  normalisation  condition  has  been  chosen  such  that  r°(z=0)  =  1 .  The 
compatibility  condition  leads  to  the  neutral  stability  relation  Ma^=8  k  (k+Bi). 
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We  can  now  define  adjoint  vectors 


(A.  10) 


belonging  to  a  set  F*  of  adjoint  boundary  conditions  to  be  defined  later  on. 
A  scalar  product  is  introduced  in  the  usual  way  (f.ex.  [14,24])  by 


{v\S^(U))  =  [  [v;.s,  +  w^  S,]dz  +  (A.ll) 

-00 

where  5,  i=l,..,5  are  the  components  of  (A. 8),  and  the  overbar  denotes  the  complex 
conjugate.  We  also  define  the  adjoint  operator  Sj  of  by 

[v*  ,SjiU))  =  {si(y*),u)  (a.12) 


for  all  U  belonging  to  E  and  V*  to  F*.  Integration  by  parts  leads  to 


siiv*) 


iD^-k^)V*  -ikp* 
(D^-k^)W*  -Dp*  +  T' 
DW*  +iT.y; 


(A.  13) 


[  (P^-k^)T*  J 

and  the  cancellation  of  the  boundary  term  gives  F*  as  the  set  of  sufficiently  derivable 
functions  satisfying 

V*  ,W*  ,DT*  ,p*  forz^-o:>  (A.14) 

W*  =  DV*  =X*+V*  =  DT*  +BiT*  +iMa^'k.X*  =0  forz=0  (A.15) 


The  resolution  of  the  adjoint  problem  Sj  (Vj )  =  0  (with  Vj  €  F  * )  gives 

ikk'^(l  -kz  -k^z^) 
z(l  -kz) 

.  =  - -  _4(1  +kz) 

*  A{2k+Bi) 

8  ^ 


(A.  16) 


where  the  normalisation  (  V-r  ,  6t(Fj)  /  =  1  has  been  adopted. 
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Finally,  the  eigenfunctions  of  the  spectral  problem  (12)  read 


' 

V- 

ikk  2e*^(l+^z) 

X 

Zgkz 

=  -4/t2(A:+B/) 

2e'‘^ 

pi 

^iz(2  ^  -  2 yte <■'**'>'''"(  ^  +  ^  ) 

Tl 

o  k^Ma 

<  > 

(A.  17) 


where  the  corresponding  eigenvalue  a  is  solution  of  the  dispersion  relation  (13).  Note  that 
the  derivation  of  (A.  17)  assumes  a+l^>0  (which  is  verified  provided  that  k<Ma/2  Bi,  as 
seen  from  eq.  (13)).  Otherwise,  when  a+k^<0,  imaginary  roots  are  obtained  from  the 
characteristic  relation,  and  lead  to  a  continuum  of  solutions  bounded  at  z-^-oo ,  but  which  do 
not  satisfy  boundary  conditions  (2). 


Appendix  3  :  Calculation  of  the  mean  temperature  profile 


For  the  roll  mode  Aj  =  a^(t)8  {k-  k^)  +  ay(t)8(k  +  k^)  (see  end  of  section  III. 3),  the  k=0 
Fourier  component  of  the  perturbation  vector  is  found  from  (A. 5)  as 

2|fl.P[/i’(z)  (A.  18) 


/t=0 


where,  according  to  (A. 6),  Uoi(z)  is  the  solution  of  the  inhomogeneous  problem 

^o(^o?)  =  (A.19) 

since  M^^g—O,  and  the  resonant  part  of  Nj  _j  is  zero,  as  ZjQ  —  O  from  (A. 3).  Here  again,  the 
subscripts  refer  to  the  mode  number  (i  for  k^,  -1  for  -^q).  It  is  computed  that  the  only  non¬ 
zero  component  of  C/o?(z)  is  Tqx{z),  solution  of 

.-/lorf,  =  z)(<7’,‘’)  (A-19) 

where  the  neutral  stability  functions  are  given  by  (A. 9).  Then,  the  solution  belonging  to  E 
is 

To?  =  (55/2  ^  nk^Bi+lkl){e^'‘^\^ -z)  -  ^)  +  (Bi+k,yky^°hH5 -2k,z)  (A.20) 


The  cubic  coefficient  Zfn 
Ziu  =  -{v: 


can  then  be  computed  from  (A. 7)  as 


))-l 


”  ,  0,  0  k^(Bi+kAH3Bi+5kA 

TiiWifTUz-  — 


2(Bi+2k^ 


(A.21) 


A  slightly  longer  but  similar  computation  leads  to  the  expression  of  the  second  cubic 

coefficient  ,  and  then  finally  to  eq.  (20).  Equation  (27)  is  also  obtained  in  this  way,  but 
by  using  eigenfunctions  (A.  17)  instead  of  the  neutral  stability  functions  (A. 9). 
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IV.  NUMERICAL  RESULTS 


TV.l.  Description  of  the  method 

In  this  section,  a  full  numerical  simulation  of  the  nonlinear  equations  is  presented.  Details 
on  the  method  are  reported  in  [1],  such  that  we  only  present  here  the  principles  of  the 
numerical  scheme.  The  two-dimensional  geometry  used  for  calculations  is  represented  on 
figure  1.  We  will  assume  insulating  and  "slippery"  lateral  walls  atx=0  and  a:=L  (where  L 
is  the  aspect  ratio  of  the  container),  which  have  already  been  used  to  model  real  rigid  walls 
[2].  This  is  equivalent  to  the  periodic  boundary  conditions  used  in  section  III.  Here,  we  will 
consider  a  bottom  plate  at  z=0,  which  is  rigid  and  heat  conducting.  The  free  surface  at  z=l 
is  assumed  undeformable  and  insulating  (Bi=0). 


The  liquid  layer  is  submitted  to  a  linear 
temperature  profile,  and  the  evolution  of 
perturbations  of  this  reference  state  is 
smdied.  The  dimensionless  equations  we 
consider  are 


—  =  Pr(AV  -  Vp  +  RaTU)  (4-2) 
Dt 

—  =  AT  +VA-  (4-3) 

Dt  " 


where  DIDt  is  the  convective  derivative  Figure  4.1  :  Geometry  of  the  smdied  system. 
DIDt  =  didt  +  (V.V)  and  V  =  {U,W),  T 

and  p  are  respectively  the  deviations  of  velocity,  temperamre  and  hydrodynamic  pressure 
around  the  mechanical  equilibrium  state  (see  section  II).  A  thermal  time  scale  h^Z/c  is  used 
as  in  the  previous  sections  (h  is  the  depth  of  the  liquid  layer). 


Note  that  the  Rayleigh  effect,  i.e.  the  buoyancy  forces  induced  by  density  variations  with 
temperamre,  have  been  included  as  the  last  term  of  eq.  (4.2).  The  dimensionless  number 
quantifying  the  importance  of  these  forces  on  the  stability  of  the  diffusive  state  is  the 
Rayleigh  number  Ra=ga/3h''/j'fc  (g  is  the  gravity  acceleration,  a  the  thermal  expansion 
coefficient,  jS  the  imposed  thermal  gradient,  v  the  kinematic  viscosity  and  k  the  thermal 
diffusivity). 

The  boundary  conditions  are  derived  from  figure  1  : 


W  =  dWIdz  =  r  =  0 

at  z=0 

(4.4) 

W  =  dT/dz  =  0 

(4.5) 

d^W/dz^  -Mad^TIdx^  = 

0  at  z=l 

U  =  dT/dx  =0  at 

x=Q,L 

(4.6) 
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where  the  usual  notations  are  used. 


The  Galerkin  method  consists  in  first  expanding  the  unknown  quantities  in  series  of  functions 
of  X  and  z  with  time-dependent  coefficients  : 

W{X,Z,t)  =  Y. 

m.n  (4.7) 

T(x,z,t)  =  Y 

m,n 


We  choose  the  basis  functions  and  in  order  to  satisfy  the  boundary  conditions  (4-6). 
Then  the  developments  are  introduced  in  equations  2  and  3,  which  are  then  projected  on  the 
basis  functions.  The  discretisation  of  the  problem  is  accomplished  by  truncating  the 
expansions  (7)  to  a  finite  number  of  terms.  The  obtained  system  is  of  the  following  form 


dA, 


dt 

~dr 


^Bjj  ) 
’By  ) 


(4.8) 


and  is  well  suited  for  computer  resolution.  The  choice  of  and  d^„  is  detailed  in  [1].  We 
just  mention  here  that  trigonometrical  functions  are  used  for  the  horizontal  dependency  :  i.e. 
w„,„=cos(mk(jx)W„(z),  B^„=^cos(mkipc)TJz),  which  satisfies  lateral  (periodic)  boundary 
conditions  at  x=0  and  x=L=ir/ko.  The  choice  of  vertical  functions  is  more  delicate,  due  to 
the  complexity  of  boundary  conditions  at  z=0  and  z=l.  The  computer  simulation  uses  a 
Runge-Kutta  method  of  order  4,  with  adaptative  stepsize  control  for  numerical  integration 
given  some  initial  conditions  B^„(0).  We  always  chosen  infinitesimal  "white  noise" 

A, JO)  =  B,J0)  =  10-^. 


TV. 2.  Finite-amplitude  regimes  of  convection  -  behaviour  at  large  Ma 

The  following  figures  represent  convective  states  obtained  for  various  values  of  the 
Marangoni  number  Ma,  the  Rayleigh  number  Ra  being  in  all  cases  equal  to  zero.  We  have 
considered  four  different  values  of  the  aspect  ratio  L  of  the  box  :  1.57,  0.785,  0.628  and 
0.524.  As  the  lateral  walls  can  be  considered  as  periodic  boundary  conditions,  these  aspect 
ratios  correspond  respectively  to  values  2,  4,  5  and  6  of  the  wavenumber  kg  of  the  convective 
field  (since  one  convective  cell  corresponds  to  one  half  of  the  period,  the  aspect  ratio  L  and 
the  wavenumber  ko  are  related  by  ko=7r/L). 

For  large  values  of  the  wavenumber  kg  (small  values  of  L),  the  convection  cells  are  more 
concentrated  near  the  interface,  and  consequently  do  not  feel  the  presence  of  the  lower  rigid 
plate.  Thus,  for  kg-^oo,  the  results  obtained  here  should  tend  to  the  results  obtained  in  the 
last  section  III,  where  the  liquid  layer  was  assumed  infinitely  deep  (short  wave  effects).  This 
provides  a  way  to  test  the  assumptions  used  in  section  III,  which  allowed  to  obtain  analytical 
results.  This  will  be  done  in  section  IV. 3. 

For  each  simulation,  the  following  data  (steady  state)  are  presented  : 

-  the  streamline  pattern  (top  left  figure,  the  direction  of  circulation  is  clockwise) 

-  the  isotherms  (top  right  figure,  on  which  the  dark  zones  indicate  cold  regions,  while  the 
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white  zones  indicate  hot  regions) 

-  the  mean  temperature  profile,  i.e.  the  temperature  field  averaged  over  one  horizontal  period 
(middle  left  figure,  where  the  abscissa  is  the  vertical  coordinate  z,  and  the  ordinate  is  the 
mean  temperature) 

-  the  free  surface  temperature  as  a  function  of  the  horizontal  coordinate  x  (middle  right) 

-  the  free  surface  horizontal  velocity  as  a  function  of  the  horizontal  coordinate  x  (bottom  left) 

For  each  figure,  numerical  values  of  typical  convective  quantities  are  also  presented  :  these 
are 

-  the  increase  of  the  mean  surface  temperature  (averaged  in  the  horizontal  direction)  due  to 
convection.  This  convective  quantity,  denoted  as  delta,  is  a  measure  of  the  increase  of  the 
heat  transport  created  by  convection,  as  remarked  in  section  III.  For  ko-><»,  it  should 
correspond  to  the  value  of  the  bulk  temperature  decrease  A  of  section  III  (see  IV. 3). 

-  the  amplitude  of  the  surface  temperature  variations  (denoted  by  delta  T),  defined  as  the 
difference  between  temperature  maxima  and  minima  on  the  free  surface.  For  ko-^oo ,  it  should 
correspond  to  the  value  of  ATs„rf  of  the  section  III. 

-  the  value  of  the  maximal  surface  velocity  Vmax,  together  with  the  horizontal  position  x  of 
this  maximum.  Again,  for  ko-»<»,  it  should  correspond  with  values  of  computed  in 
section  III. 

Also  given  on  each  figure  are  the  Marangoni  number  Ma,  the  value  of  the  wavenumber  ko, 
the  number  Nmax  of  vertical  functions,  and  the  number  Imax  of  horizontal  trigonometric 
functions  in  the  Galerkin  expansions  (4.7). 

For  each  value  of  kg,  the  Marangoni  number  is  increased  step  by  step  until  the  point  where 
numerical  divergences  occur.  These  divergences,  due  to  the  insufficient  number  of  trial 
functions  in  the  Galerkin  expansion,  manifest  first  in  the  curve  of  the  free  surface  velocity 
(bottom  left).  This  will  also  be  apparent  on  the  curves  of  the  free  surface  maximal  velocity 
presented  in  the  section  IV. 3.  However,  some  of  these  slightly  diverged  simulations  are 
presented  here,  and  it  will  be  seen  in  the  next  section  that  the  values  obtained  for  other 
convective  quantities  may  still  be  considered  as  satisfactory. 
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IV.3.  Comparison  of  numerical  with  analytical  results 

In  order  to  compare  results  of  sections  HI  and  IV,  we  first  have  to  remark  that  the  length 
scales  used  to  put  the  problem  under  dimensionless  form  are  different.  In  section  III,  we 
choosed  an  arbitrary  length  d  to  scale  the  problem,  while  in  section  IV,  the  depth  h  of  the 
layer  was  choosed  for  this  purpose.  As  the  analytical  results  of  section  III  were  obtained  for 
a  value  k=l  of  the  dimensionless  wavenumber  (which  means  the  dimensional  value  of  this 
wavenumber  is  k*=  1/d),  the  corresponding  value  of  the  dimensionless  wavenumber  in  section 
IV  is  k’.h=h/d=ko  (value  given  for  each  numerical  simulation). 

The  velocities  were  scaled  by  x/h  in  section  IV,  and  by  K/d  in  section  III.  The  dimensional 
value  V’  of  the  velocity  being  independent  of  the  choice  of  units,  we  must  have 
V*=V‘'.K/d=V'’.ic/h  (where  V*  and  V”  stand  for  dimensionless  values  of  the  velocity  obtained 
respectively  in  section  III  and  IV).  Thus,  the  transformation  rule  between  dimensionless 
velocities  is  V‘‘=V''/ko.  Similarly,  as  the  Marangoni  number  Ma  is  proportional  to  the  square 
of  the  unit  of  length,  we  have  the  equivalence  Ma‘'=Ma''/ko2.  Finally,  Ae  unit  of  temperature 
being  /3d  (section  III)  or  jSh  (section  IV),  the  temperature  equivalence  is  T*=T''.ko. 


Figure  4.2  :  The  increase  of  the  mean  surface  temperature  due  to  convection  (in  units  of 
i8d=i3/k*),  i.e.  the  bulk  temperamre  decrease  as  a  function  of  the  Marangoni  number.  Plain 
curves  ;  numerical  results  (the  value  of  kg  is  indicated  on  each  curve).  Thick  dashed  curve: 
eq.  3.28,  thin  dashed  curve  :  the  limit  of  eq.  3.28  for  Ma-^oo,  i.e.  8/3—2.66.  A  good 
agreement  is  obtained  for  ko-»><» ,  as  in  this  limit,  the  convection  cells  are  localised  near  the 
interface  (see  previous  figures),  and  the  layer  can  effectively  be  assumed  infinitely  deep. 
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Figure  4.2  presents  the  values  of  the  increase  of  the  mean  surface  temperature ’due  to 
convection,  expressed  in  units  of  section  III  (thus,  values  obtained  in  section  IE  are 
unchanged,  while  numerical  results  of  section  IV  have  to  be  multiplied  by  ko  :  this  gives 
delta. ko),  as  a  function  of  the  Marangoni  number  Ma=-aT]8dV/iK=-ffTi8//xKk*2.  The  thick  full 
curves  correspond  to  numerical  results  of  section  IV,  corresponding  to  the  values  2,4,5,  and 
6  of  the  wavenumber.  The  thick  dashed  line  corresponds  to  the  analytical  formula  (3.28)  of 
section  III  for  the  bulk  temperature  decrease  (which  is  equivalent  to  the  free  surface 
temperature  increase),  and  the  thin  dashed  horizontal  line  to  the  limit  of  this  expression  for 
Ma-^oo,  i.e.  8/3  =  2.66  (as  given  by  eq.  3.40).  It  is  seen  that  all  the  numerical  curves  seem 
to  present  a  saturation  of  the  surface  temperature  increase  when  Ma-»oo,  although 
calculations  would  be  necessary  at  higher  Marangoni  numbers,  but  would  require  a  much 
larger  computational  effort.  Furthermore,  the  numerical  curves  approach  the  analytical  ones 
when  ko  is  increased,  as  expected  from  the  discussion  at  the  beginning  of  section  IV.2.  In 
fact,  the  analytical  results  obtained  in  section  III  appear  to  be  quite  reliable  for  computing 
this  convective  quantity  in  the  limit  ko-^oo ,  which  is  a  very  satisfactory  result  taking  into 
account  the  drastic  assumptions  underlying  their  derivation  (mean-field  approximation, 
limitation  of  amplitude  equations  to  cubic  terms  and  negligence  of  harmonics  of  the 
fundamental  wavenumber). 


Figure  4.3  :  The  maximal  surface  velocity  (in  units  /c/d=/ck*)  as  a  function  of  the  Marangoni 
number  Ma=-CTT/3/)Lt/ck*^.  Plain  curves  :  numerical  results  (the  value  of  ko  is  indicated  on  each 
curve).  Thick  dashed  curve  :  analytical  (using  3.27),  thin  dashed  curve  :  asymptotic 
behaviour  for  Ma-»-oo,  i.e.  3.67  Ma‘'^  (eq.  3.41).  Again,  a  good  agreement  is  obtained  for 
ko-»oo,  apart  from  numerical  inaccuracies  causing  a  divergence  of  the  curves  above  Ma=40. 
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The  behaviour  of  other  convective  quantities  with  Ma  is  represented  on  figs  4.3  (maximal 
surface  velocity)  and  4.4  (amplitude  of  surface  temperature  variations).  It  is  seen  that  the 
agreement  is  also  quite  satisfactory,  already  for  values  of  ko=4  and  5. 


Figure  4.4  :  The  amplitude  of  the  surface  temperature  variations  (in  units  i8d=|8/k*)  as  a 
function  of  the  Marangoni  number  Ma=-aj^/iJLKk*^.  Plain  curves  ;  numerical  results  (the 
value  of  ko  is  indicated  on  each  curve).  Thick  dashed  curve  :  analytical  result  (using  eq. 
3.27),  thin  dashed  curve  :  asymptotic  behaviour  for  Ma-»oo,  i.e.  14.66  (eq.  3.42). 
Again,  a  good  agreement  is  obtained  for  ko-*oo,  already  for  values  of  ko=4  and  5. 
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rv.4.  Conclusions  of  section  IV 


In  this  section,  a  full  numerical  integration  of  the  equations  governing  Marangoni-Benard 
convection  in  finite  idealised  boxes  has  been  presented.  Although  the  Galerkin  formulation 
includes  the  effect  of  buoyancy,  the  value  of  the  Rayleigh  number  as  been  considered  as  zero 
(microgravity,  or  layers  with  small  depth).  In  this  first  analysis  of  the  problem,  the  free 
surface  Biot  number  has  been  set  to  zero.  Extension  to  non-zero  Biot  number  (e.g.  due  to 
evaporation)  could  be  achieved  in  the  future. 

The  lateral  walls  of  the  layer  have  been  assumed  slippery  (which  means  that  the  velocity 
tangential  to  the  wall  is  not  zero,  as  would  be  the  case  for  real  rigid  walls),  and  insulating. 
This  kind  of  walls,  which  have  often  been  used  in  previous  studies  of  the  effect  of  lateral 
walls  on  nonlinear  convection,  are  known  to  simplify  the  mathematical  analysis  of  the 
problem.  In  particular,  they  are  equivalent  to  periodic  boundary  conditions,  thus  allowing 
straightforward  comparison  with  analytical  results  of  section  III. 

Numerical  simulations  have  been  presented  for  values  of  the  wavenumber  Iq,  equal  to  2,4,5, 
and  6  (corresponding  respectively  to  values  L=1.57,  0.785,  0.628  and  0.524  of  the  aspect 
ratio  L=l/h  of  the  box,  where  1  and  h  are  respectively  the  dimensional  width  and  height  of 
the  box).  For  each  value  of  ko,  the  Marangoni  number  is  increased  step  by  step.  For  each 
step,  the  following  results  are  presented  for  the  steady  convective  state  reached  after  transient 
have  been  damped  out :  the  streamline  pattern,  the  isotherms,  the  mean  temperature  profile, 
the  free  surface  temperamre,  and  the  free  surface  horizontal  velocity.  Numerical  values  of 
typical  convective  quantities  are  also  presented  :  these  are  the  increase  of  the  mean  surface 
temperature  due  to  convection  (which  is  equivalent  to  the  bulk  temperature  decrease  A  of 
section  III,  and  measures  the  enhancement  of  heat  transfer  due  to  convection),  the  amplitude 
of  the  surface  temperature  variations,  the  value  of  the  maximal  surface  velocity,  together 
with  the  horizontal  position  of  this  maximum. 

The  results  show  that  the  enhancement  of  heat  transfer  due  to  Marangoni  convection  presents 
a  saturation  behaviour  when  the  Marangoni  number  is  increased.  The  free  surface  velocities 
grow  continuously  with  Ma,  and  the  amplitude  of  the  surface  temperature  variations 
(difference  between  hot  and  cold  spots  on  the  free  surface)  first  grows  and  then  decreases 
when  Ma  is  increased.  For  ko-»oo ,  it  is  seen  that  numerical  results  compare  very  well  with 
analytical  results  obtained  via  the  amplitude  equations  analysis  of  section  III,  despite  the 
drastic  assumptions  making  its  realisation  possible.  This  is  a  proof  of  the  validity  of  these 
approximations  for  studying  Marangoni-Benard  convection  far  from  the  instability  threshold. 
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V.  GENERAL  CONCLUSIONS 

In  the  first  part  of  this  work,  we  first  computed  the  rate  of  evaporation  of  a  liquid  layer  lying 
on  a  heated  rigid  plate  at  constant  temperature  Tb,  the  free  surface  of  which  being  in  contact 
with  a  gas  phase  at  constant  pressure  Pg.  Evaporation  is  modelled  taking  into  account 
deviations  from  equilibrium  (existence  of  a  chemical  potential  difference  on  both  sides  of  the 
free  surface,  revealing  a  deviation  of  the  state  of  the  interface  from  the  Clausius-Clapeyron 
coexistence  curve),  which  may  occur  for  example  in  the  presence  of  small  quantities  of 
impurities,  imperfect  gases,  ...  resulting  in  low  values  of  the  accommodation  coefficient 
appearing  in  &e  Hertz-Knudsen  law  (equivalently,  low  values  of  the  phenomenological 
coefficient  K  appearing  in  our  formulation  of  the  problem).  This  first  analysis  resulted  in  the 
determination  of  a  particular  solution  (the  basic  state)  of  the  full  system  of  conservation 
equations  (mass,  momentum  and  energy)  describing  the  problem.  A  graphical  determination 
of  the  evaporative  mass  flux  J  was  proposed,  together  with  an  anal5d:ical  formula,  valid 
provided  the  energy  flux  (proportional  to  the  mass  flux)  is  not  too  large  (of  the  order  of  some 
Watts/cm^  for  a  1mm  depth  of  water). 

In  a  second  stage,  we  determined  the  conditions  under  which  such  a  basic  state  is  stable 
against  hydrodynamic  fluctuations,  allowing  the  surface  tension  of  the  liquid  to  depend  on 
temperature.  Indeed,  evaporation  at  the  free  surface  creates  a  temperature  gradient  inside  the 
liquid  layer  near  the  interface,  such  that  this  situation  can  become  unstable  above  a  certain 
value  of  this  gradient.  In  the  limit  of  small  liquid  Peclet  numbers  (relatively  small 
evaporative  mass  fluxes),  the  temperature  profile  is  linear  in  the  liquid  phase.  If  furthermore 
the  Crispation  number  is  small  (and  the  Galileo  number  is  large),  the  surface  deformation  can 
be  neglected,  and  the  stability  problem  is  shown  to  reduce  to  the  classical  Pearson's 
problem.  The  critical  value  of  the  Marangoni  number  Ma  (proportional  to  the  value  of  the 
mass  flux  and  to  the  surface  tension  variation  with  temperature)  is  determined  as  a  function 
of  the  Biot  number,  characterising  the  effective  heat  transfer  through  the  free  surface,  and 
function  of  the  evaporation  parameters.  In  particular,  it  can  be  computed  that  this  number 
can  be  as  large  as  10^  for  typical  liquids  as  water.  However,  this  value  can  be  reduced  by 
one  or  two  order  of  magnitudes  by  the  presence  of  unavoidable  impurities  on  the  free 
surface,  reducing  the  values  of  the  accommodation  coefficient  (in  the  Hertz-Knudsen 
relation),  and  accentuating  the  deviation  from  the  equilibrium  limit.  The  critical  Marangoni 
number  is  found  to  increase  with  the  Biot  number,  thus  revealing  a  stabilising  influence  of 
evaporation. 

A  nonlinear  analysis  of  the  convection  problem  was  then  undertaken  in  section  III,  in  order 
to  determine  the  enhancement  of  heat  transfer  created  by  convection,  when  the  threshold  of 
stability  is  exceeded.  This  analysis  is  based  on  the  assumptions  that  the  dynamics  of 
convection  can  be  described  by  the  interactions  of  the  unstable  eigenmodes  of  the  problem, 
that  the  layer  can  be  assumed  infinitely  deep  (this  is  not  contradictory  with  the  assumption 
of  a  small  Peclet  number,  provided  that  attention  is  restricted  to  short  wave  effects,  localised 
near  the  free  surface),  and  that  the  Prandtl  number  of  the  liquid  can  be  assumed  infinite.  It 
was  found  that  thermocapillary  convection  results  in  a  decrease  of  the  bulk  temperature  (or 
equivalently  in  an  increase  of  the  free  surface  temperature),  such  that  the  injected  thermal 
energy  is  transported  under  a  smaller  temperature  difference,  due  to  the  additional  heat 
transport  created  by  convection.  The  value  of  this  bulk  temperature  decrease  was  found  to 
increase  with  the  Marangoni  number,  and  shown  to  present  a  saturation  behaviour  when  the 
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Marangoni  number  tends  to  infinity.  Analytical  formulas  were  obtained,  as  a  function  of  the 
Biot  number.  Other  convective  quantities  have  also  been  computed  (the  maximal  surface 
velocity,  behaving  as  Ma*'^  for  Ma->oo ,  and  the  amplitude  of  surface  temperature  variations, 
decreasing  as  Ma‘^'^  in  the  same  limit).  Results  were  also  obtained  concerning  the  nonlinear 
competition  between  unstable  modes  :  in  particular,  the  wavelength  selection  problem  has 
been  investigated  on  the  basis  of  a  numerical  transient  integration  of  the  system  of  amplitude 
equations.  It  was  shown  that  in  a  first  stage,  convection  always  develops  from  infinitesimal 
perturbations  under  the  form  of  convective  structures  having  horizontal  wavelengths 
determined  by  the  fastest  growing  perturbations  (wavenumbers  proportional  to  Ma*'^  for  large 
Ma).  When  these  perturbations  become  sufficiently  large,  convective  cells  begin  to  coalesce, 
and  to  form  larger  and  larger  wavelength  structures.  This  process  becomes  slower  and 
slower,  due  to  the  increasing  inertia  of  the  fluid  to  set  in  motion,  and  finally  stabilises  when 
the  size  of  convection  cells  becomes  comparable  to  the  size  of  the  vessel  containing  the 
experiment.  The  bulk  temperature  decrease  is  found  to  increase  continuously  during  this 
process,  thus  indicating  the  general  tendency  of  natural  convection  to  increase  the  heat 
transfer.  Other  results  were  obtained  about  the  competition  between  roll  and  hexagonal 
structures.  It  was  shown  that  the  hexagonal  structures  are  generally  preferred  for  moderate 
Marangoni  numbers,  while  rolls  stabilise  for  very  large  Marangoni  numbers.  This  confirms 
earlier  results  qualitatively,  even  if  the  method  (and  the  quantitative  results)  are  significantly 
different. 

Finally,  a  full  two-dimensional  numerical  simulation  of  the  problem  was  performed,  by  the 
use  of  a  Galerkin  method.  This  allowed  to  incorporate  the  effect  of  a  finite  depth  of  liquid 
(lower  rigid  plate).  Results  were  obtained  for  various  aspect  ratios  of  the  convective  cells  by 
varying  the  aspect  ratio  of  the  simulation  domain.  In  particular,  a  very  satisfactory 
(quantitative)  agreement  has  been  obtained  with  analytical  results  obtained  via  the  amplitude 
equations  analysis,  when  the  aspect  ratio  (width/depth)  tends  to  zero,  in  which  case  the 
convective  perturbations  are  effectively  located  near  the  interface,  and  do  not  perceive  the 
presence  of  an  horizontal  rigid  bottom  plate.  This  rather  good  concordance  is  a  proof  of  the 
legitimacy  of  the  above-mentioned  assumptions,  which  should  thus  be  helpful  for  further 
smdies  of  Marangoni-Benard  convection  heat  transfer  characteristics. 
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